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Numerical analysis of influence of ship hull form 
modification on ship resistance and propulsion 
characteristics 


Part II 
Influence of hull form modification 
on wake current behind the ship 


Tomasz Abramowski, Ph. D. 
Tadeusz Szelangiewicz, Prof. 
West Pomeranian University of Technology, Szczecin 


ABSTRACT 


After signing ship building contract shipyard s design office orders performance of ship resistance and 
propulsion model tests aimed at, apart from resistance measurements, also determination of ship speed, 
propeller rotational speed and propulsion engine power for the designed ship, as well as improvement of 
its hull form, if necessary. Range of ship hull modifications is practically very limited due to cost and time 
reasons. Hence numerical methods, mainly CFD ones are more and more often used for such tests. In this 
paper consisted of three parts, are presented results of numerical calculations of hull resistance, wake 
and efficiency of propeller operating in non-homogenous velocity field, performed for research on 18 hull 
versions of B573 ship designed and built by Szczecin Nowa Shipyard. 


Keywords: ship hull geometry, numerical (computational) fluid dynamics, resistance, wake, propeller efficiency 


NUMERICAL COMPUTATIONS of model tests [1] are given for comparison. The average value 

OF WAKE CURRENT of wake fraction determined numerically within the range of 

values of the relative propeller radius r/R reaching from 0.25 

The second effect of ship hull form modification, apart to 1.21 is equal to 0.48, while that obtained from the model 
from change of resistance, is change of wake current which tests amounts to 0.52. 

to a large extent influences performance of screw propeller 

(its thrust, torque and efficiency) as well as overall propulsive 

efficiency of ship. 
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T 
Model basin | 
[deg «e [deg] 
Fig. 14. Circumferential distribution Fig. 16. Circumferential distribution 
of axial wake fraction component for r/R — 0.202 of axial wake fraction component for r/R — 1.0 
In Fig. 17 stream lines for the inital ship hull form are also 
shown. 
ọ [deg | 
Fig. 15. Circumferential distribution 
Fig. 17. Distribution of stream lines on ship hull model 


of axial wake fraction component for r/R = 0.60 


Tab. 5. Results of numerical calculations of the axial wake fraction component W. [-] 
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Tab.6. Results of model tests of the axial wake fraction component W. [-] 


RESULTS OF ANALYSIS OF WAKE 
CURRENT FOR MODIFIED VERSIONS OF 
SHIP HULL 


Numerical calculations of wake current were performed 
for the same modified ship hull versions as in the case of the 
resistance investigations (Tab. 2). 

The ship hull model speed V = 1.492 [m/s]. 

The collected results of numerical calculations in the form 
of the average wake fraction according to Eq. (2), are given 
in Tab. 7. 

The influence of particular modified geometrical parameters 
of ship hull model on the average value of wake fraction is 
presented in Figs. 18 + 20, and, in Fig. 21 and 22 are given 
circumferential distributions of axial and tangential wake 
fraction components for r/R = 0,6. 


Wake fraction component Wx 


Hull block coefficient Cp 
Fig. 18. Wake fraction component W, calculated in function of C, 


Number of 
variant 


Tab. 7. Average values of axial wake fraction 
for modified hull versions of B 573 ship 


Modified 
parameter 


Average value of axial 
wake fraction 


51% 
Manual 
modification 
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Wake fraction component Wy 
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Hull prismatic coefficient C, 


Fig. 19. Wake fraction component W calculated in function of C, 
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component Wx 
Ux) ON 


Wake fraction 
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Longitudinal location of centre of buoyancy LCB 
[in percent of waterline length measured from FP] 


Fig. 20. Wake fraction component W, calculated in function of LCB 
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Fig. 21. Circumferential distribution of field of the axial wake fraction 
component for r/R = 0.6, calculated for all the analyzed variants 
(18 in number) 
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Fig. 22. Circumferential distribution of field of the tangential velocity 
component for r/R — 0.6, calculated for all the analyzed variants 
(18 in number) 
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In Figs. 23 + 29 are exemplified distributions of 
tangential and axial wake fractions as well as velocity 
vectors on the propeller disk plane for selected variants of 
B 573 ship hull modification. The complete set of numerical 
calculation results is contained in the report on the research 
project [2]. 
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Fig. 23. Distributions of tangential velocities and axial wake fraction 
as well as velocity vectors in propeller disk plane 
for variant 1 of B 573 ship hull modification 
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Fig. 24. Distributions of tangential velocities and axial wake fraction Fig. 25. Distributions of tangential velocities and axial wake fraction 
as well as velocity vectors on propeller disk plane for variant 6 as well as velocity vectors on propeller disk plane for variant 7 
of B 573 ship hull modification of B 573 ship hull modification 
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Fig. 26. Distributions of tangential velocities and axial wake fraction 
as well as velocity vectors on propeller disk plane for variant 11 
of B 573 ship hull modification 
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Fig. 27. Distributions of tangential velocities and axial wake fraction 
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of B 573 ship hull modification 


Vr [m/s] 


ü 20 40 60 80 100 120 140 160 180 
Angle [deg] 
Distribution of the tangential velocity V, 


Wy 


0 20 40 60 80 100 120 140 160 180 
Angle [deg] 
Distribution of the axial wake fraction component W, 


1.67e+00 
1.59e+00 
1.51e«00 
1.42e+00 
1.34e+00 
1.26e-00 
1.17e+00 
1.09e+00 
1.01e+00 
9.22e-01 
8.39e-01 
7.55e-01 
6.72e-01 
5.88e-01 
5.05e-01 
4.21e-01 
3.38e-01 
2.54e-01 
1.71e-01 
8.75e-02 Ls 


4.04e-03 


Wake current velocity vectors 


Fig. 28. Distributions of tangential velocities and axial wake fraction 
as well as velocity vectors on propeller disk plane for variant 17 
of B 573 ship hull modification 
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ABSTRACT 


After signing ship building contract shipyards design office orders performance of ship resistance and 
propulsion model tests aimed at, apart from resistance measurements, also determination of ship speed, 
propeller rotational speed and propulsion engine power for the designed ship, as well as improvement of 
its hull form, if necessary. Range of ship hull modifications is practically very limited due to cost and time 
reasons. Hence numerical methods, mainly CFD ones are more and more often used for such tests. In this 
paper consisted of three parts, are presented results of numerical calculations of hull resistance, wake 
and efficiency of propeller operating in non-homogenous velocity field, performed for research on 18 hull 
versions of B573 ship designed and built by Szczecin Nowa Shipyard. 


Keywords: ship hull geometry, numerical (computational) fluid dynamics (CFD), 
resistance, wake current, propeller efficiency 


SHIP PROPULSION EFFICIENCY 


In design practice the overall propulsive effciency is presented 
in the form of the product of efficiency components [1]: 


N = Nao lars Ne (3a) 
or: 
N = Neo res" (3b) 
where: 
Nar — ship hull ,,efficiency” determined under the 
assumption that the open-water -propeller thrust T 
is equal to the behind-the-hull propeller thrust T 
1 — ship hull ,,efficiency” determined under the 


assumption that the torque Q, delivered to the open- 
water propeller is equal to the torque Q,, delivered 
to the behind-the-hull propeller 

Tlo — open-water propeller efficiency 

rotative ,efficiency" determined under the same 

assumption as for hull ,,efficiency” 

Ns — shaftline efficiency 

Ng — reduction gear efficiency. 


Mep Tro — 


The two components: n, and n,, are mechanical efficiencies, 
not associated with water flow around ship’s hull, propeller 
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and rudder. Therefore they do not affect ship’s hull-propeller 

- rudder collaboration. 

The first three components, being of hydrodynamic 
character, are decisive of overall propulsive efficiency, and 
basin model tests to be performed have to provide ship 
designer with information necessary a.o. to determine their 
values. 

However worth paying attention is that : 

e for determining n, and n, are used averaged values 
obtained from measurements taken during model basin 
tests (the wake current velocity field measured during 
model tests is easily converted to ship scale, that not always 
brings good results) 

* the presented efficiency components do not contain 
other factors affecting the overall propulsive efficiency 
(influence of some of them can be determined by 
conducting multiple model basin tests e.g. for various 
versions of ship hull stern part) 

* in the above presented formulas is contained the open- 
water-propeller efficiency n, determined for the screw 
propeller operating in homogenous water velocity field 
without presence of ship hull, whereas the behind-the-hull 
screw propeller operates in a very non-homogenous velocity 
field, and the coefficients n,,, n, do not reflect essence of 
the problem, 


* also, in the above presented formulas effects of rudder, 
especially streamline one, which may increase efficiency 
of screw propeller, are lacking 

e the above presented formulas do not contain information 
on mutual location of screw propeller and blade rudder 
relative to ship stern, and on its impact on propulsive 
efficiency. 


The above described way of the determining of propulsive 
efficiency contains certain ambiguities (e.g. those associated 
with wake-fraction determining) hence in the publications [2, 3] 
it is postulated to introduce the coefficients connected with the 
behind-the-hull propeller efficiency: 


_ Qo 
so On 


_T 
er ni T 
called the torque modifier and thrust modifier, respectively. 


(4a) 


(4b) 


Regardless of that which approach would be applied to 
determining the ship propulsive efficiency, it is rational to 
determine thrust and torque of the behind- the-hull screw 
propeller as well as its efficiency in non-homogenous water 
velocity field for mean working point equivalent to mean 
service speed of ship (the mean statistical ship speed was 
discussed in [4]). 

To determine real propulsive efficiency of ship or also 
efficiency of propeller operating behind the hull in non- 
homogenous velocity field, is possible by using CFD techniques. 
Hence in this part of the paper are presented results of numerical 
computation of efficiency of screw propeller operating in non- 
homogenous water velocity field, performed for 18 modified 
hull versions of B 573 ship. For the computations use is made of 
the results of numerical calculation of wake currents, presented 
in the part II of the paper. 


AOR. 


0.9R 


NUMERICAL CALCULATION OF 
HYDRODYNAMIC CHARACTERISTICS 
OF SCREW PROPELLER 


Before the calculations of efficiency of a screw propeller 
operating in non-homogenous velocity field, hydrodynamic 
characteristics of a store screw propeller used to propulsion 
model tests of B 573 ship [5], were calculated (the screw 
geometry is presented in Fig. 30). 

The calculations of the store screw propeller were 
performed in the scale 1:1 by modelling the flow with the use 
of RANS equations and RNG k -epsilon model was applied 
to turbulence phenomena. Boundary layer was modeled by 
means of standard flow approximating functions, where the 
parameter y* = 50. 

To discretize domains a numerical non-structural network 
of prismatic elements within boundary layer (Fig. 31), which 
was further adjusted to polyhedral elements (Fig. 32). was 
applied. The non-structural network was elaborated by using 
Gambit system; and the conversion was performed by using 
Fluent system. All screw propeller blades were modeled, that 
demands on one hand a greater computer capacity, but on the 
other hand makes non-stationary calculations of screw propeller 
with taking into account non-design conditions, possible. 

Calculation time step was equal to 0.00012 s, that was 
equivalent to about 1? of one screw-propeller revolution at 
its rotational speed of 17 rps. Two complete screw-propeller 
revolutions were covered by the calculations, hence the total 
time of the calculated flow was equal to 0.0864 s. 


Results of numerical calculations: 

— pressure distribution on screw propeller blades is presented 
in Fig. 33 and 34 

— streamlines — in Fig. 35 

— total velocity profiles — in Fig. 36 

— hydrodynamic characteristics — in Fig. 37 + 39 
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Fig. 30. Geometry of the store screw propeller No. P 355 [5] 
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Fig. 33. Pressure distribution on driving face of screw propeller blades, 
the advance ratio J = 0.4 


Fig. 34. Pressure distribution on drag surface of screw propeller blades, the advance 
ratio J = 0.4 
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Fig. 36. Total velocity profiles 
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Fig. 37. Calculation results of the thrust coefficient KT 
compared with model test results [5] 
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Fig. 38. Calculation results of the torque coeficient KQ 
compared with model test results [5] 
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Fig. 39. Calculation results of the screw propeller efficiency eta 
compared with model test results [5] 


MEAN EFFICIENCY OF FINAL SCREW 
PROPELLER IN NON-HOMOGENOUS 
WAKE CURRENT VELOCITY FIELD 


The modification of B 573 ship’s hull form was aimed at 
investigations of influence of the introduced hull form changes 
on the ship’s resistance, wake current and screw propeller 
efficiency. The obtained distributions of wake current were 
used for determining the mean efficiency of the final screw 
propeller for B 573 ship. 

The numerical calculations were performed for the final 
screw propeller of B 573 ship [6] in the same way as for the store 
screw propeller, however in the case of the final screw propeller, 
the water velocity field (in which axial and circumferential 
components were taken into account), calculated for each of 
the modified version of the ship’s hull form, was introduced. 
The final results of the analyses are presented in Fig. 40 in the 
form of the mean efficiency of screw propeller operating in 
non-homogenous wake current field. In the figure the mean 
efficiency values are shown for particular versions (variants) 
of the modified ship hull forms in compliance with Tab. 1. The 
case 0 is calculated for the screw propeller operating in the non- 
homogenous flow velocity field measured during model tests 
of the B 573 ship hull (of non-modified form). The continuous 
line (marked red) stands for the maximum efficiency of the 
final open-water screw propeller operating in homogenous 
flow velocity field. 


0 12 3 4557 83 9 10 11 12 13 
Variants of modified hull forms of B 573 ship 


14 15 16 17 18 19 


Mean efficiency of the final screw propeller for B 573 


Fig. 40. Mean efficiency of B 573 ship 5 screw propeller operating in non- 
homogenous wake current field, where: y, — maximum efficiency of the final 
screw propeller in homogenous flow velocity field, variant 0 — original hull 
form (without modification) in wake current velocity field measured during 

model tests, variant 18 — symmetrical hull form with manually modified 
stern part 


FINAL CONCLUSIONS 


- The tests performed in advance each ofthe main calculations 
demostrated a good conformity of CFD calculations to 
model tests. Certain observed differences may result from 
various causes (some of them were discussed in the Part I of 
this work). The differences in hydrodynamic characteristics 
of the store screw propeller could also result from some, 
out of all necessary to elaborate a 3D model and numerical 
network, geometrical quantities lacking in [5]. Therefore it 
is not certain whether the computational model of the store 
screw propeller has represented exactly its real features. In 
the case of the final screw propeller there was no available 
results of model tests or other calculations to make 
verification of the CFD computation results possible. 

- Numerical calculations of the mean efficiency of screw 
propeller for 18^ hull version of B 573 ship (Fig. 40) 
demonstrated that small modifications can have either 
favourable or unfavourable effect onto efficiency. Especially 
favourable effect in the form of elevated efficiency was 
obtained by using manual modification of only stern part 
of the ship's hull (18" variant, Fig. 40). It means that by 
suitable modelling hull form in its stern part only it is 
possible to make flow velocity field before the propeller 
more favourable from the point of view of efficiency. 

- Mean values of screw propeller efficiency for wake currents 
corresponding to particular versions of ship hull were 
calculated for the same final screw propeller. Investigations 
on searching for optimum screw propeller geometry for 
each of the ship’s hull versions, are under way. 

- Numerical methods based on the CFD provide even greater 
possibilities, e.g. investigations on propulsive efficiency of 
the whole system composed of ship hull, screw propeller 
and streamline rudder, that would make performing a more 
advanced optimization of ship propulsion system, possible. 
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ABSTRACT 


The computer system for the complete design of the contra-rotating propellers presented in this article has 

several common blocks and procedures with the systems for design of single propellers and tandem co- 

rotating propellers, presented in detail in the Polish Maritime Research No.1 and No.4 of the Volume 16, 

2009. In this article only the blocks and procedures developed specially for the contra-rotating propellers 

are described. The system is based on the lifting line and lifting surface models and on the Computational 

Fluid Mechanics methods. The comparative analysis of the contra-rotating propellers and the tandem co- 
rotating propellers is included. 


Keywords: ship propellers, contra-rotating propellers, design methods, computational fluid dynamics 


INTRODUCTION 


Many years of experience and practice have shown that 
the classical single propellers in most cases fulfil the highest 
requirements regarding efficiency, cavitation characteristics, 
level of unsteady hydrodynamic shaft forces, pressure 
pulsations generated on the hull and noise/vibration. However, 
in realistic applications certain areas may be found in which 
other types of screw propellers are better. Consequently, the 
contra-rotating propellers (CRP) may be advantageous in 
some specific conditions. For example, the requirement of 
minimum slipstream rotation is achievable for single shaft 
installations only with contra-rotating propellers (a very 
important requirement for underwater vehicles). 

The design of contra-rotating propellers is based on the 
same requirements and assumptions which are employed in the 
design of single propellers and tandem co-rotating propellers 
[4, 6]. In general, the design requirements of contra-rotating 
are similar to those of tandem propellers, but they differ in the 
following characteristic features: 

a) in the design program for contra-rotating propellers the 
determination of the aft propeller diameter is required 

b) detailed information about the velocity field induced by 
both propellers of the set is required not only in the design 
program, but also in the program for analysis of the contra- 
rotating propellers operation in the non-uniform wake. 


The long-term practical experience and the analysis of the 
propeller-induced velocity fields shows that: 

- the diameter of the aft propeller must be correlated with 
the diameter of the slipstream behind the forward propeller 
in such a way that the tip vortices shed from the forward 
propeller pass outside the tips of the aft propeller blades 
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- the number of blades of the forward and aft propellers of 
the contra-rotating set should be different 

- the rate of rotation of both propellers of the contra-rotating 
set is recommended to be different. 


Similarly as in the case of design of other propeller types the 
computer design system for contra-rotating propellers includes 
three interacting programs (blocks of procedures): 

1) program for determination of the design velocity field for 
both propellers of the contra-rotating set 

2) program for propeller design 

3) program for analysis of the contra-rotating propellers 
operation in the design velocity field, including the effects 
of mutual interaction between both propellers. 


In comparison to the system for design of single propellers 
the system for design of contra-rotating includes the following 
special elements (similarly as the system for tandem 
propellers): 

a) procedures for the independent graphical presentation of 
both propellers of the set 

b) procedures for determination of the design velocity field 
for both propellers, taking into account their mutual 
interaction 

c) procedures for modification of the velocity field for the 
analysis of the contra-rotating propeller operation in the 
non-uniform inflow velocity field 

d) procedures for determination of the induced pressure 
pulsations for each propeller separately and for the entire 
set of contra-rotating propellers 

e) procedures for determination of the unsteady hydrodynamic 
forces and moments for each propeller separately and for 
the entire set of contra-rotating propellers. 


iosa fato ala ques aee rd 


Fig. 1. The block diagram of the computer system for the complete design of the contra-rotating propellers 
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The appropriate co-operation between all programs and 
procedures, supplemented with the above listed elements, 
ensures effective and complete design of the contra-rotating 
propellers. The system is organized in such a way that the 
designer controls the entire design process from the computer 
screen and the data transfer between the respective elements 
of the system is fully automatic. 

Similarly as in the case of the previously presented systems 
[4, 5, 6] an extensive use has been made of computer graphics 
for control and correction of the input data, for control of the 


cum 


intermediate results, for modification of the geometry of both 
propellers at consecutive stages of design and for presentation 
of the final results. 

The block diagram of the system for design of the contra- 
rotating propellers, shown in Fig. 1, is similar to that for 
design of tandem propellers [6]. The basic block diagram of 
the propeller design process is extensively supplemented with 
blocks for the design of contra-rotating propellers (cf. Fig. 2) 
and for the analysis of contra-rotating propellers operation in 
the non-uniform inflow velocity field (cf. Fig. 3). 


PropantToDes 


(INPUT DATA) 


DETERMINATION OF VELOCITY 
INDUCED BY FORWARD PROPELLER 
ON AFT PROPELLER 


P2 4 ^ EPSS0.005 


| DETERMINATION OF GEOMETRY AND 
PRESSURE DISTRIBUTION ON BLADES 
OF AFT AND FORWARD PROPELLER OF 


Fig. 2. The block diagram of the design procedure of the contra-rotating propellers 
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Fig. 3. The block diagram of the analysis procedure of the contra-rotating propellers 


PRESENTATION OF THE SELECTED 
BLOCKS OF THE NEW COMPUTER 
SYSTEM 


The computer system for design of the contra-rotating 
propellers has several blocks identical to those incorporated 
in the systems for design of single or tandem co-rotating 
propellers. Only the blocks concerning the specific problems 
of contra-rotating propellers are described in detail below. 


The input data 


The input data include all information required for initiation 
of the variant of design calculation selected out of four 
available options, similarly as for other types of propellers 
described in earlier publications [4, 5, 6]. The input data may 
be introduced either in the form of the previously prepared 
input data file or directly from the computer terminal in an 
interactive mode. 

The input magnitude which is typical both for the tandem 
co-rotating and contra-rotating propellers is the axial distance 
between the generator lines of the forward and aft propeller. The 
features distinguishing the contra-rotating propellers from the 
tandem co-rotating propellers is the difference in diameters of 
the forward and aft propellers and the recommended difference 
in number of blades and rate of rotation. The diameter of the 


aft propeller must be smaller than the local diameter of the 
slipstream shed from the forward propeller. This requirement is 
based on two reasons: firstly, there is a marked jump in velocity 
at the slipstream boundary (see Figs 9 and 10), which cannot 
be accommodated by the blade geometry (it would require 
a corresponding jump in the blade pitch) and secondly, there 
is danger of cavitation erosion caused by the tip vortices shed 
from the forward propeller. 

The program enables graphical control and correction of 
the input data, especially the radial distributions of geometrical 
parameters such as blade outline, maximum thickness of the 
blade section profiles, blade skewback, inflow velocity in the 
ship wake etc. Fig. 4 shows an example of correction of the 
radial distribution of the circumferentially averaged inflow 
velocity in the ship wake. 


The design program 


The algorithm of the design program for contra-rotating 
propellers differs in the following details from the analogical 
programs for single or tandem co-rotating propellers [4, 6]: 

a) the calculations are performed only for the given value of 
the total thrust of the CRP set 

b) the division of the total thrust between the forward and 
aft propeller may be given in the input data or may be 
determined in the design program 
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Fig. 4. The control display of the given radial distribution of the circumferentially averaged velocity in the ship wake, 
showing the correction of the erroneously set value 


c) the diameter of the aft propeller may be given in the input 
data or may be determined in the design program (the 
determined value may still be corrected by the designer 
during calculations) 

d) the determination of the velocity field induced by both 
propellers of the CRP set is necessary (this is done 
automatically in the design program). 


The design calculations are performed in the same way for 
every selected version ofthe design task (similarly as in the case 
of classical single propellers). The results of design calculations 
are available in the form ofthe appropriate numerical files and 
drawings on the computer screen, which may be saved and, or 
printed. An example of such a drawing is shown in Fig. 5. This 
drawing may be viewed on the screen from different angles and 
then printed in the most suitable configuration. 


è 
Pi 


Fig. 5. The contra-rotating propellers as represented in the design program 


As it has been mentioned above, the determination of 
the velocity fields induced by the bound and trailing vortex 
systems representing both propellers of the CRP set is an 
indispensable element ofthe design program. It may be seen in 
the block diagram shown in Fig. 2 that the design of CRP is an 
Iterative process based on the mutual hydrodynamic interaction 
between the both propellers of the CRP set. The results of these 
calculations are presented in the format similar to that of the 
tandem co-rotating propellers [6]. 
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The induced velocity field 


In the propeller design program based on the vortex theory 
the systems of bound and trailing vortices are determined 
in a simplified form, which is sufficient for calculation of 
the velocities induced on the blades of the propellers, but it 
may be insufficient for calculation of the velocities induced 
at certain distance in front and behind the propeller blades 
(i.e. at locations of the respective propellers of the CRP or 
tandem co-rotating set of propellers). 

The trailing vortex system of any propeller undergoes 
contraction, deformation and concentration (rolling-up of 
vortices with simultaneous viscous dissipation of vorticity) 
[7, 8, 9, 10]. As a result, tip vortices are formed behind every 
blade and a strong hub vortex is formed behind the propeller 
hub. The velocity field induced by such a system of concentrated 
vortices is strongly three-dimensional, being a function of all 
three co-ordinates (for example x, ®, r in the cylindrical 
system). In the design task of CRP set we are interested in the 
circumferentially averaged values of the induced velocities, 
while in the case of tandem co-rotating propellers the fully 
local induced velocity field must be known. The examples of 
the induced velocity field behind a propeller in the selected 
cross-section x — const ofthe slipstream are shown in Figs 6, 7 
and 8. The velocities induced by the aft propeller at the location 
of the forward propeller have much smaller values and more 
uniform spatial distribution. 

In case of the CRP set the velocity field induced by the 
forward propeller, combined with the inflow velocity field 
of the hull wake is an important design parameter. All three 
velocity components of that induced velocity field are averaged 
circumferentially (on the basis ofthe flow volumetric intensity) 
for a required number of radii. The radial distribution of these 
averaged velocity components shows a marked discontinuity at 
the local tip vortex radius. At the distances behind the propeller 
larger than the propeller radius both velocity components 
(axial and tangential) are equal zero outside the slipstream 
and they reach quite high values inside the slipstream. Such 
a sudden jump of the velocity cannot be accommodated by the 
appropriate geometry of the aft propeller blades. Hence the 
requirement to reduce the aft propeller diameter to the value 
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Fig. 6. The axial component of the velocity 
induced behind a three-bladed propeller 
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Fig. 7. The tangential component of the velocity 
induced behind a three-bladed propeller 
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Fig. 8. The radial component of the velocity 
induced behind a three-bladed propeller 


which ensures that the tips of the blades are already inside the 
slipstream shed from the forward propeller. 

The volumetric averages of the velocities induced in front 
of the aft propeller are much smaller, their radial distribution 
is smooth and it does not exhibit any discontinuities (the 
volumetric average of the tangential induced velocity in front 
of the propeller is theoretically equal zero). 


The program for analysis of the contra-rotating 
propellers operation in the non-uniform 
velocity field 


The computer program UNCA for the analysis of the ship 
propeller operation in the non-uniform velocity field is a very 


important element of the design process of all propeller types. The 
central part of the algorithm of this program is the determination 
of the size and extent of various forms of cavitation on the blades 
of the propeller operating in the circumferentially non-uniform 
velocity field. The original theoretical model combines the 
unsteady vortex lifting surface model with the model of unsteady 
sheet cavitation bubble. The detailed description of the algorithm 
of UNCA may be found in [18, 19]. 

In the case of a CRP set the modification of the non-uniform 
velocity field of the ship hull wake due to the velocity field 
induced by the propeller other than this currently analyzed is 
a very important part of the computation algorithm. In order to 
fulfil this task an additional procedure has been incorporated 
in the program This procedure calculates the induced velocity 
field both in front and behind the propeller. 

The appropriately modified velocity field enables the detailed 
analysis of the operation of both forward and aft propeller by 
means of the program UNCA. In comparison with other propeller 
types the analysis of the CRP operation in the non-uniform 
inflow is particularly complicated and time-consuming. The aft 
propeller of the set operates in the non-uniform velocity field 
being the result of superposition of the hull wake and the velocity 
field induced by the forward propeller (see Figs. 9, 10 and 11). 
This analysis is additionally complicated by the fact that for each 
angular position of the forward propeller all selected angular 
positions of the aft propeller must be analyzed. The total number 
of the analyzed angular blade positions depends also on the 
number of blades of both propellers of the CRP set, but as a rule 
this number exceeds one thousand (for example for Zf= 5 and Za 
= 4 the total number of the analyzed angular blade positions is 
n = 35*36 = 1260). For each of these blade positions the 
appropriate numerical files containing results are generated, 
which undergo further re-calculations, but they also may be 
presented on the computer screen in the numerical or graphical 
form or directly printed. These results include, among other 
data, the calculated cavitation phenomena on the aft propeller, 
shown in Fig. 12. 

The results of calculations of the unsteady hydrodynamic 
forces and moments on the aft propeller for the angular 
position of the forward propeller equal to AF = 0 are presented 
below (Figs. 13, 14 and 15). Similar diagrams may be 
presented for every mutual position of the forward and aft 
propeller. Independently, analogical results are available for 
the forward propeller. 
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Fig. 9. The velocity field behind the ship hull 
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Fig. 10. The velocity induced by the forward propeller at the location of the 
aft propeller at the angular position of the forward propeller AF = 0 
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Fig. 11. The resultant velocity field at the location of the aft propeller for the 
angular position of the forward propeller AF = 0 
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Fig. 12. The example of computed cavitation phenomena 
on the contra-rotating propellers 


20 POLISH MARITIME RESEARCH, No 1/2010 


40 so $00 700 
Angle [deg] 


Fig. 13. The computed hydrodynamic force components FX, FY, FZ on the 
aft propeller for the angular position of the forward propeller AF = 0 
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Fig. 14. The computed hydrodynamic moment components MX, MY, MZ on 
the aft propeller for the angular position of the forward propeller AF = 0 


For every angular blade position of the forward propeller 
of CRP the maximum and minimum values of the calculated 
pressure pulsations generated on the hull by the aft propeller 
and of the unsteady shaft forces are extracted. These data are 
presented in the numerical and graphical form (see Figs. 17 
and 21). 

Both propellers of the CRP set are analyzed by means of 
the program UNCA from the following points of view: 

a) the intensity and extent of the different forms of cavitation 
in different angular positions of both propellers in the 
combined non-uniform velocity field including the hull 
wake and the propeller-induced velocities, 

b) the level of pressure pulsations generated by the CRP set 
and by the forward and aft propeller separately on the hull 
surface or in the surrounding space, 

c) the level of the fluctuating hydrodynamic shaft forces and 
moments for the forward and aft propeller separately. 


After such an analysis the designer may return to the 
CRP design calculation, modifying the selected details of the 
geometry and hydrodynamic parameters of both propellers of 
the set: 

a) the radial distribution of blade skewback 

b) the radial distribution of the blade section chord lengths 

c) the radial distribution of the maximum blade section 
thickness 

d) the type of blade section thickness distribution 

e) the type of blade section mean line 

f) the radial distribution of hydrodynamic loading of the 
blades 

g) the number of blades of both propellers 

h) the division of CRP total thrust between the forward and 
aft propeller 

1) the diameter of one or both propellers of CRP. 


The above described analysis of CRP is performed only 
for the design condition, because any change of the ship 
speed or the propellers rates of rotation changes the operating 
conditions of both propellers — consequently the fields of 
induced velocities of both propellers change in the qualitative 
and quantitative sense. 

The numerical results of the above analysis are stored 
in three separate files, referring respectively to the forward 
propeller, aft propeller and to the entire CRP set. 

The method of graphical presentation of the results is similar 
to that applied to the classical single propellers [4], with the 
modification that in case of CRP it may refer separately to the 
forward or aft propeller or to the entire CRP set. The examples 
of such presentation are shown in the next section. 


COMPARATIVE ANALYSIS OF CRP AND 
TANDEM CO-ROTATING PROPELLERS 


In the preceding publications [4,6] the results of design 
calculations of the single propeller and of the tandem co- 
rotating propellers. The case of the same large fast cargo vessel 
has been used in the design of the CRP set described below. The 
design ship speed is V = 25.3 [knots] and the required thrust 
of the propulsor is T = 3250 [kN] The design rate of rotation 
of the engine and of the forward propeller is 100 [rpm]. In the 
case of CRP set the rate of rotation of the aft propeller may 
be different and it may be optimized from the point of view 
of propulsive efficiency, cavitation, pressure pulsations on the 
hull and unsteady hydrodynamic shaft forces. In the process 
of such an optimization the best performance of CRP was 
achieved for the number of blades of the forward propeller 
Zf = 5, number of blades of the aft propeller Za = 3 and for 
the rate of rotation of the aft propeller equal to 115 [rpm]. The 
optimum diameter of the aft propeller was determined on the 
basis of the appropriate semi-empirical formulae concerning 
the slipstream behind the forward propeller. The calculated 
main characteristics of the forward and aft propellers of the 
CRP set are given in Tab. | below. 


Tab. 1. Main parameters of the designed CRP set 


Forward 
propeller 


7.80 
0.59 
30060 
23343 


Characteristic Aft propeller 


Propeller diameter [m] 


Expanded blade area ratio [-] 


Shaft power [kW] 
Mass of the blades [kg] 


Static moment of inertia 
[kgm**2] 


511944 


The comparison of the above results for CRP with those 
obtained for other propeller types [4, 6] shows that the propulsive 
efficiency of the CRP set is about 5 percentage points higher than 
that ofa single propeller and about 2.5 percentage points higher 
than the tandem co-rotating propellers. However, the CRP set is 
worse than the other analyzed propellers as far as the pressure 
pulsations on the hull (probably also acoustic pressures) and the 
unsteady hydrodynamic shaft forces are concerned. 

Figs. 15 and 16 below present the calculation results for the 
CRP set, concerning the harmonic amplitudes of the unsteady 
hydrodynamic shaft forces and moments, shown separately for 
the forward and aft propeller. Presentation of the results is done 
separately for each of the propellers of the CRP set because 
their different numbers of blades and different rates of rotation 
prohibit presentation of the results on one common diagram of 
harmonic amplitudes. 


Moreover, the calculations of the harmonic amplitudes of the 
unsteady hydrodynamic shaft forces are performed separately 
for each angular position of the forward and aft propellers. 
The numerical files containing the results or their graphical 
representations may be inspected. However, the analysis of 
about 1000 diagrams or files is not easy. In order to facilitate 
this analysis, for each angular position of one of the propellers 
only the maximum and minimum values of the respective 
components of the unsteady shaft forces and moments are 
extracted. These values are appropriately added together and 
presented in the form of a relatively short numerical file or 
as one diagram. Example of such a diagram is presented in 
Fig. 17. This diagram includes also the differences between the 
maximum and minimum values because these differences are 
decisive parameters in some Classification Society Rules. 


Fig. 15. The computed harmonic amplitudes of the hydrodynamic force 
and moment components for the forward propeller of the CRP set 


| FY FZ Mx MY MZ 
Fig. 16. The computed harmonic amplitudes of the hydrodynamic force 


and moment components for the aft propeller of the CRP set 
for the angular position AF = 0 of the forward propeller 


-4 000.00 
Fig. 17. Maximum and minimum values, together with their differences, 
for the hydrodynamic shaft forces [kN] and moments [kNm] 
on the entire CRP set 
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For the sake of comparison the corresponding results of 
calculations for the tandem co-rotating set, designed for the 
same ship and operating in the same non-uniform velocity field 
[6], are presented in the analogical format in Fig. 18. 


Fig. 18. Maximum and minimum values, together with their differences, 
for the hydrodynamic shaft forces [RN] and moments [kNm] on the set of 
tandem co-rotating propellers 


A similar presentation of the calculation results may be 
done for the pressure pulsations generated on the hull (or 
acoustic pressures in the surrounding space). In Figs. 19 
and 20 the harmonic amplitudes of the pressure pulsations 
generated in the selected points on the hull by each propeller 
of the CRP set separately are presented. Fig. 21 shows the 
maximum and minimum values of these pressure pulsations, 
together with their differences. For the sake of comparison, 
the corresponding values calculated for the tandem co-rotating 
propellers, designed for the same ship and operating in the same 
non-uniform velocity field are shown in Fig. 22. 


pressure 
588 


Harmonic amplitudes of 


Point number 


Fig. 19. Harmonic amplitudes of the pressure pulsations generated 
on the hull by the forward propeller of the CRP set 


Similarly as in the case of other propeller types, the program 
determines the presence and extent of the various cavitation 
forms, which are described as numerical files or corresponding 
illustrations. In Figs. 23, 24 and 25 the examples of graphical 
presentation of the computed cavitation phenomena are shown. 

The cavitation pictures may be viewed independently for 
each of the propellers (Figs. 23 and 24) or in combination 
(Fig. 25). Apart from the static pictures the program enables 
viewing the CRP set in motion, exhibiting the calculated 
unsteady cavitation phenomena. The complete revolution of 
both propellers of the CRP set may be presented, showing 
the time-dependent variation of the cavitation phenomena 
on the propeller blades in the non-uniform velocity field (in 
this presentation the forward or aft propeller is alternatively 


22 POLISH MARITIME RESEARCH, No 1/2010 


immobilized and the other performs one full revolution). This 
presentation may be also recorded in the format of a film. 


Harmonic amplitudes of pressure p 


Point number 


Fig. 20. Harmonic amplitudes of the pressure pulsations generated 
on the hull by the aft propeller of the CRP set 


Fig. 21. Maximum and minimum values, together with their differences, 
of the pressure pulsations [kPa] induced on the hull by the entire CRP set 


Fig. 22. Maximum and minimum values, together with their differences, 
of the pressure pulsations [kPa] induced on the hull by the entire set 
of tandem co-rotating propellers 


Tip vortex 
cavitation 
Fig. 23. Example of the computed cavitation phenomena 

on the forward propeller of the CRP set 


‘Intermittent 
cavitation 


Tip vortex 
cavitation 


Fig. 24. Example of the computed cavitation phenomena 
on the aft propeller of the CRP set 


Bubble 
cavitation 


Tip vortex 
cavitation 


Fig. 25. Example of the computed cavitation phenomena on both propellers 
of the CRP set for their selected position in the non-uniform inflow 


It may be concluded from the analysis of the above 
examples, that from the point of view of propulsive efficiency 
the CRP set is better than the tandem co-rotating propellers and 
than the classical single propeller. This is not a general rule, the 
design example presented in [4,5,6] is selected in such a way, 
that it may be demonstrated that in certain areas of propeller 
operation the CRP or tandem co-rotating propellers may be 
more efficient than the classical single propeller. 

The problem of cavitation performance, unsteady 
hydrodynamic shaft forces and pressure pulsations is much 
more involved and it requires a thorough analysis of the 
results of comparative calculations. The computer systems 
presented in the current line of publications [4, 5, 6] facilitate 
the optimum selection of the appropriate propeller type for 
every application. 


CONCLUSIONS 


e The computer system presented above fully facilitates the 
process of design of the contra-rotating propellers. This 
system incorporates all elements necessary in the correct 
design process of CRP. The computational process itself is 
relatively fast (although much slower than in the case of 
other, simpler propellers), and the graphical subroutines 
enable efficient analysis of the consecutively designed 
CRP sets, directing the whole process towards an optimum 
solution of the design task. 


e The relatively short computation time allows for calculation 
of a large number of variants of the designed CRP set. For 
example the following parameters may be modified in the 
course of design without leaving the system: 

- number and geometry of the blades of both propellers 
of the CRP set (including the geometry of the blade 
sections) 

- diameters of both propellers of the CRP set 

- ship speed and rates of rotation of both propellers of the 
CRP set (equal or different) 

- division of the total thrust between the forward and aft 
propeller of the CRP set 

- radial distribution of the hydrodynamic loading of the 
blades 

- axial distance between both propellers of the CRP set. 


e The above listed modifications enable the optimization of 
the CRP set from the point of view of efficiency, cavitation 
performance, induced pressure pulsations and unsteady 
hydrodynamic shaft forces and moments. The widely 
employed graphical presentation of the results greatly 
facilitates such an optimization analysis. 


e The above presented design example, together with the 
results of calculations of the ducted propeller included in [5] 
and tandem co-rotating propellers included in [6], lead to the 
conclusion that it is worthwhile to include in the design 
analysis different propeller types. The four propeller design 
systems presented in the current line of publications [4, 5, 6], 
which accept very similar input data files (relatively the 
most onerous part of the design process), enable very fast 
and efficient comparative analysis, leading to the selection 
and consecutive design of the optimum propeller type in 
every practical application. 
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Influence analysis of changes of design 
parameters of passenger-car ferries on their 
selected sea-keeping qualities 


Tomasz Cepowski, Assoc. Prof. 
Szczecin Maritime University 


ABSTRACT 


The main scientific aim of this research was to elaborate design guidelines which could 

make it possible to improve sea-keeping qualities of passenger-car ferries. The searched- 

for design guidelines were prepared in the form of regression functions as well as artificial 

neural networks on the basis of the results obtained from calculations with the use of 

numerical methods based on the theory of planar flow around a body. The guidelines make 

it possible to predict ship roll, sea-sickness index, lateral and vertical accelerations on the 
basis of quantities available in the preliminary stage of ship design. 


Keywords: passenger-car ferry, design guidelines, preliminary design stage, sea-keeping qualities, rolling, 
sea- sickness (motion-sickness) index, vertical accelerations, lateral accelerations, ship main dimensions 


INTRODUCTION 


Many dividing criteria of sea-going ships can be found 
because of variety of floating units, e.g. according to their 
functions, size, propulsion, architectural form, design problems, 
spatial subdivision etc. For most types of ships sea-keeping 
qualities are not an important design limitation. Exceptions 
from the rule are certain types of ships which can be divided 
into the following groups: 

1. group of ships which have to fulfil their mission independent 
of weather conditions 

2. group of ships characterized by design features which 
increase their susceptibility to weather conditions 

3. group of ships intended for the carrying of passengers 

4. group of transport ships for which sea-keeping qualities are 
only one of the limitations imposed on them. 


The ships belonging to the groups have to be characterized 
by various sea-keeping qualities which can be modeled on 
different design levels (preliminary design — technical design). 
For each of the above mentioned ship groups important sea- 
keeping qualities resulting from the functions to be fulfilled by 
them, can be formulated. 

In this work only the group of passenger ships has been 
selected out of all the groups. The group can be further 
divided into: 

1. ships intended for the carrying of passengers only: 

- excursion ships (cruisers) 

- passenger ships 

- sanatorium and hospital ships 

- yachts 


2. ships intended for the carrying of passengers and general 
cargo: 
- passenger-cargo ships 
- passenger-cargo ferries: 
passenger-car ferries 
passenger-train ferries 
- ropax ships 
- passenger-container carriers. 


The main design problems of ships of the group are first of 
all: their internal capacity, speed, and necessity of providing 
safety and comfort for passangers. 
Among the first subgroup of the ships a. o. excursion 
ships (cruisers) are numbered. And, with a view of their large 
dimensions, to provide them with good sea-keeping qualities 
is not a problem. However for the remaining ships of the 
group, especially passenger ships intended for carrying cargo, 
design solutions are searched for to improve their sea-keeping 
Tunes such as these aimed at: 
the decreasing of roll motion which can cause cargo to shift 
and in consequence to impair ship stability 
* the decreasing of accelerations which can cause lashing of 
fixed cargo to break (or cargo to shift) 

* the moderating of accelerations which can generate sea- 
sickness 

* the avoiding of parametric resonance of motion. 


The passenger-car ship is a typical representative of the 


ships’ group regarding the problems associated with the 
providing-for of appropriate sea-keeping qualities. 
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PROBLEM OF MODELING SEA-KEEPING 
QUALITIES OF PASSENGER-CAR FERRIES 
DURING PRELIMINARY DESIGN STAGE 


In ship design process design solutions which satisfy both 
economic criteria and technical limitations, are searched for. 
In the case of passenger-car ferry the economic criteria are 
consisted of many requirements put by shipowner, to which 
a.o. internal hull capacity or service speed can belong. And, the 
technical limitations are formed from a.o. a group of factors 
affecting ship safety, i.e. first of all appropriate ship stability, 
unsinkability and hull strength. 

Ship designers are often and often required to design some 
of the passenger-car ferries so as to additionally provide them, 
apart from general design merits, with immunity to weather 
conditions, i.e. good sea — keeping qualities. 

Hull form and dimensions of a ferry are to a great extent 
decisive of its sae-keeping qualities. Therefore they should be 
modelled already in the preliminary design stage. An important 
feature of preliminary ship design is that exact form of ship 
hull is represented by means of the so-called main dimensions 
and global coefficients which characterize hull form, e.g. hull 
block, midship-section and vertical prismatic coefficients.The 
scarce set of information makes using the known methods of 
determination ship motions in waves, based on the classical 
linear or non-linear ship motion theory, not possible; hence in 
present it is important problem how to take fully into account the 
entire range of sea-keeping qualities already in the preliminary 
ship design stage. And, another problem is that the selection 
of improper values of ship main dimensions and hull form 
coefficients irreversibly worthens sea-keeping qualities (and 
other design merits as well) since any change of any dimensions 
of a ship after its building is economically unjustified. 

Therefore this work has been aimed at the elaborating of 
design guidelines focused on improvement of sea-keeping 
qualities of sea-going passenger-car ferry already in the 
preliminary stage of its design. 


METHOD OF INVESTIGATIONS 


To reach the aim of this research the method consisted of 

the following phases, has been selected: 

1. elaboration of a list of ferry-ships covering a broad range 
of their forms and dimensions 

2. elaboration of operational scenarios important from the 
point of view of ferry operation and significantly affecting 
certain technical features of such ships 

3. simplification of physical model by: 
a. replacement of real weather conditions by statistical 

ones 


b. replacement of real waves by a standard wave energy 
spectral density function 

c. replacement of instantaneous ship response values by 
conventional statistical ones. 

4. elaboration of symbolic mathematical model which 
describes the ferry design parameters X, X,, ... X , as well 
as the sea-keeping qualities K important from the point of 
view of an assumed ferry operational scenario 

5. use of the exact numerical methods to calculate values of 
the parameter K 

6. selection of approximating functions for a set of discrete 
values of K-parameter, and elaboration of design guidelines 
with taking into account the parameter in question 

assessment of the elaborated design guidelines 

determination of a range of application and limitation of 
the elaborated formulas. 

The detailed research alg 


DA 


orithm is presented in Fig. 1. 


List of variants 
of hull forms and dimensions 
of passenger-car ferries 


Conventional 
operational 
scenario 


Ferry motion 
and wave 
parameters 


Calculation of K-paramter 
model values by using 
the exact methods 


Analysis of the results 
and determination 
of the parameters X,, X,. ... X, 
which significantly 
influence the parameter K 


Elaboration of design guidelines 


Fig. 1. Algorithm of the applied research method, where: X, X, ... X, — ship 
design parameters, K — a parameter which describes ship sea-keeping 
qualities 


LIST OF VARIANTS OF HULL FORMS 
AND DIMENSIONS OF PASSENGER-CAR 
FERRIES 


While elaborating the list of variants of hull forms of 
passenger-car ferries the recommendations contained in the 
report [4] were taken into account. In this research the list 
of 3072 variants was prepared by using a combination of the 
following items: 


Tab. 1. Variants of passenger-car ferry hull form 


CB(Lpp) 
[-] 


Variant 


CWL(Lpp) XB/Lpp XF/Lpp 


[%] 


[-] [%] 


CB - block coefficient of underwater ship hull part, CB(Lpp) — block coefficient of underwater ship hull part related to ship length between perpendiculars, 
CM - midship section coeffcient, CB(L) — cylindrical coefficient of underwater ship hull part, CB(V) — vertical prismatic coefficient of underwater ship hull 
part, CWL — waterplane coefficient, CWL(Lpp) — waterplane coefficient related to ship length between perpendiculars, XF — distance between geometrical 
centre of waterplane and aft perpendicular, XB — distance between centre of buoyancy and aft perpendicular, Lpp — ship length between perpendiculars 
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a) the ranges of ferry dimensions: e B/d=3,3.5,4,4.5 
e LBd = 19000, 28 000, 37 000, 46 000 mì (where: b) the set of eight hull form variants (Fig. 2), represented by 


L — waterline length, B — waterline breadth, d — ship the global hull form coefficients contained in Tab. 1 
draught) c) the values of the initial transverse metacentric height GM, 


* L/B=5.8, 6.6, 7.4, 82 ranging from 0.4 m to 1.4 m every 0.2 m. 


Variant 7 Variant 8 
Fig. 2. Hull forms of passenger-car ferries, used in the research in question 
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OPERATIONAL SCENARIOS 


In the research it was assumed that the design guidelines 
will be elaborated for ferries operating in conventional 
service conditions. It was the conditions important from 
the point of view of its operational merits, in which a ferry 
could operate. Such approach has been based on the target 
assignment concept recommended by IMO Maritime Safety 
Committee for the elaborating of new assessment standards 
of ship stability [5]. 

Ship sea-keeping qualities are significantly infuenced by 
a.o. wave parameters, ship motion parameters (ship speed 
and wave encounter angle) as well as hydrodynamic ship hull 
parameters dependent on a form of its underwater part and 
weight distribution. For instance in Fig. 3 is presented impact 
ofthe characteristic wave frequency on statistically significant 
roll amplitudes of a ship, depending on varying values of initial 
transverse metacentric height of the ship moving with steady 
parameters of its motion. 


r 
EAA 
5 10 15 20 25 
T [s] 
— GM =0.4m —GM-06m === GM=0.8m 
GM = 1.0m GM=1.2m err GM=14m 


Fig. 3. Significant ship roll amplitudes in function of the characteristic wave 
frequency T, and: GM = var, Hs = 3 m, B = 60° (where: 0° — head wave, 
180° — aft wave), v = 0 m/s 


As results from Fig. 3, the characteristic wave frequency 
together with the initial transverse metacentric height greatly 
influence ship rolling motion. From the point of view of design 
problems it is difficult to determine waving conditions in which 
a given ship operates, and especially troublesome to determine 
the characteristic wave frequency. 

Therefore only such values of the characteristic wave 
frequency which generate maximum values of the considered 
sea-keeping qualities, were taken into account in order to 
eliminate influence of their impact on ship’s behaviour. 

In Fig. 4 for instance, is presented the relation between the 
initial transverse metacentric height and significant ship roll 
amplitudes for: 

e a given value of the characteristic wave frequency (series 1) 
* varying values of the characteristic wave frequency, in 
which significant ship roll amplitudes reach their maxima 

(series 2). 


As results from Fig. 4, the series 2 show another trend and 
higher values of statistical ship rolling than the series 1. Hence 
it results that the series 2 show influence ofthe initial transverse 
metacentric height on rolling motion of ship sailing in waves 
of the most unfavourable characteristic frequency. 

With a view of the above mentioned aspects the following 
conventional operational scenarios were taken into account: 
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613 [°] 


0 0.5 1 1.5 


GM [m] 
—1—= series 1 (T = 17.7 s) —1—— series 2 (T = var) 
Fig. 4. Significant ship roll amplitudes in function of the initial transverse 
metacentric height GM: series 1 — for the characteristic wave frequency 
T = 17s, series 2 — for the characteristic wave frequency T for which 
maximum ship roll values occur, Hs = 3 m, B = 60°, v = 0 m/s 


1. the ferry has lost its propulsion (ship speed v=0 m/s) and 
found itself in waving conditions in which dangerous values 
of its rolling motion occur: 

a. the ferry is obliquely situated relative to the incoming 
wave, that results in occurrence of maximum values of 
its rolling motion 

b. the wave is of the characteristic frequency which 
produces maximum values of ship rolling motion 

c. the significant wave height Hs is contained within the 
range of 1 +3 m. 

2. the ferry sails with the service speed V = 12.5 m/s in waving 
conditions which intensify occurrence of sea-sickness 
among crew members and passengers: 

a. the ferry is obliquely situated respective to the incoming 
wave, that results in occurrence of maximum values of 
the sea-sickness index 

b. the wave is of the characteristic frequency which 
produces maximum values of the sea-sickness index 

c. the significant wave height Hs is contained within the 
range of 1 +3 m. 

3. the ferry sails with the service speed V — 12.5 m/s in 
waving conditions which generate large vertical and lateral 
accelerations in the point P (Fig. 5), which could cause 
trailers placed on car deck to shift [8]: 

a. the ferry is situated at the angle against incoming wave, 
most unfavourable from the point of view of vertical 
and lateral accelerations 

b. the wave is of the characteristic frequency which 
produces maximum values of vertical and lateral 
accelerations 

c. the significant wave height Hs is contained within the 
range of 1 +3 m. 


Fig. 5. The coordinates of the point P located on the car deck, (0.75Lpp, 
0.45B, d), for which vertical and lateral accelerations were calculated, 
where: Lpp — ship length between perpendiculars, B — ship breadth, d — ship 
draught 


NUMERICAL CALCULATIONS 


For every variant of the passenger-car ferry sailing in the 
assumed operational conditions the following sea-keeping 
qualities (the parameter K, acc. Fig. 1), were calculated; 

1. the maximum values of the significant statistical ship roll 

amplitudes, x 
2. the maximum significant values of the sea (motion) - 

sickness index, MSI „ (acc. [6]) 

3. the maximum values of the significant amplitudes of vertical 

accelerations, a and lateral ones, a in the point P, 

acc. Fig. 3. 


vmax? tmax? 


For calculations of ship motion in statistical waves the 
JONSWAP wave spectrum was assumed to be a wave energy 
spectral density function. The calculations were performed by 
using the numerical methods of the SEAWAY software, based 
on the theory of planar flow around a body. The method given 
in [3] was used for calculations of hydrodynamic coefficients. 
The accuracy tests of the SEAWAY software presented in [2], 
show a high accuracy of calculations carried out with the use 
of the program in question. 


SELECTION OF APPROXIMATING 
FUNCTIONS 


This research phase was aimed at elaboration of simplified 
analytical relations suitable for predicting selected sea-keeping 
qualities on the basis of ship design parameters. 


To preliminarily identify the parameters which significantly 
impact the assumed sea-keeping qualities, the subject-matter 
literature, analysis of parameter influence on particular sea- 
keeping qualities, as well as sensitiveness analysis was used 
by applying artificial neural networks. Next, the parameters 
characterized by the following features, were selected: 

1. suitably high variability 

2. strong correlation relative to dependent variable 

3. weak correletion relative to the remaining independent 
variables of considered model 

4. availability in preliminary ship design stage. 


The selection ofthe approximating function was performed 
by verifying first the simplest linear functions. In the case when 
the verifying hypothesis was erroneous, non-linear functions 
(power, exponential, rational ones, etc), as well as functions 
represented by artificial neural networks which reveal the 
highest accuracy of approximation and the lowest dispersion, 
were searched for. 

The approximation accuracy was assessed by using Pearson 
correlation coeffcient, and the dispersion — by means of 
standard deviation of error. 


DESIGN GUIDELINES 


The performed research resulted in the elaboration of 
the design guidelines presented in the form of the following 
relations (1), (2), (3), (4) and (5): 


1 
z e ((0-072B-1.378,11.373CB-6.532,GM-0.4}A;-A) x Å, 3.3488] +0.059 
sia I 
Dune i 0.0218 (1) 
5695 B 
© /3max = Hs : 1.6221 + - : = 0.0997 ` Ta) (2) 
VGM 
3 
MSI... = 97287997. E A 
Fw 
H 
AV. /3max 736.57: È (4) 
VFw 
1 
1+ e ([0-9071WL-929,0.072B-1.378,0.017-XF-0.926,32.787.CM-31.269.GM-0.4}:A,-As) x A, 1.107 |+ 0.066 
Alana T H; . 6) 
0.059 
where: XF -— distance between geometric centre of waterplane and 
Dia = significant roll amplitude (maximum value) [°] its aft end [m] 


] ,4 Motion Sickness (sea-sickness) 
Index (maximum value) [96] 


aV fam Significant amplitude of vertical accelerations 
(its maximum value) [°] 

at sma © Significant amplitude of lateral accelerations 
(its maximum value) [°] 

B — ship breadth [m] 


CB  - block coefficient of underwater part of ship hull [-] 
CM  - midship section coeffcient [-] 

GM - initial transverse metacentric height [m] 

Fw — waterplane surface area [m] 

LWL - waterline length [m] 


H, — significant wave height [m] 
A, — matrix of weight values: 
| -8.2648 -4.4049 2.6160 0.9067 
0.7385 . 0.6195 1.2891 -0.4663 
| -13.5773. -7.7457 -0.0425 — 1.7976 
A, -— vector of threshold values: 
[-8.2648 -4.4049 2.6160 0.9067] 
A, — column vector of weight values: 


[-3.2109 5.8162 -1.7407 5.4637] 


POLISH MARITIME RESEARCH, No 1/2010 29 


A, — matrix of weight values: 


1.244 0.400 1.205 -0.763 -1.714 
-3.021 -0.137 3.079 -4.673 -0.040 
-2.501 -0.516 0.898 -2.687 1.653 
-0.416 7.030 -0.718 0.707 5.162 

| -3.699 -0.436 -3.545 4.176 -0.364 | 


A, — vector of threshold values: 

[1.610 4.976 0.700 -2.982 4.313] 
A. — column vector of weight values: 

[2.141 2.494 1.493 1.221 -3.467] 


The relations (1) and (5) are in the form of MLP artificial 
neural networks whose structures are shown in Fig. 6. 


a) 


Fig. 6. Structures of MLP artificial neural networks used to approximate: 
a) significant roll amplitudes (maximum values), by inputs: CB, GM, B; 
b) significant amplitudes of lateral accelerations (maximum values), 
by inputs: LWL, B, XE CM, GM 


ASSESSMENT OF THE ELABORATED 
DESIGN GUIDELINES 


As results from Tab. 2, the relations (1), (2), (3), (4) and (5) 
are characterized by relatively high accuracy of approximation 
and low dispersion. The diagrams of Fig. 7 through 20 show 
comparison of accuracy ofthe formulas relative to their model 
values, as well as influence of the independent variables on 
approximation of the dependent variables. As results from 
Tab. 2 and the above mentioned diagrams all the relations are 
characterized by high correlation and low dispersion. Especially 
the relations (1) and (4) are characterized by high accuracy of 
approximation. 

18 


dwt /3max [deg] 


0 2 4 6 8 10 12 14 16 18 
ò 1/3max [deg] 


Fig. 7. Comparison of approximation accuracy relative to model values: 
model value of roll motion; 9 ,,;,,,,- roll motion value approximated 
by using Eq. (1) 
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Tab.2. Selected regression statistics 


Pearson 
regression 
coefficient R 


Standard 
deviation of | 0.48° 


Ow 3max [deg] 


0 2 4 6 8 10 12 14 16 
Pimax [deg] 
Fig. 8. Comparison of approximation accuracy relative to model values: 
model value of roll motion; 9 ,;,,,,- oll motion value approximated 
by using Eq. (2) 


Qurama È 


—— artificial neural network - Eq. (1) 


— linear regression - Eq. (2) 
+ model values 


B [m] 
Fig. 9. Comparison of approximation accuracy relative to model values: 
B = var, CB =0.61, GM = 0.8 m, H, 3m 


O usmax [deg] 


— artificial neural network - Hn (1) 


—- linear regression - Eq. (2) 
+ model values 


0.55 0.56 O57 0.58 0.59 06 O61 0.62 063 064 0.65 
CB [-] 
Fig. 10. Comparison of approximation accuracy relative to model values: 
CB = var, B = 21.6 m GM = 1.0 m, H,= 3m 


— artificial neural network - Eq. (1) 


—- linear regression - Eq. (2) 
+ model values 


0 0.2 04 0.6 0.8 1 1.2 14 1.6 
GM [m] 
Fig. 11. Comparison of approximation accuracy relative to model values: 
GM = var, B = 29.0 m, CB = 0.63, H= 3m 
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—artificial neural network - Eq. (1) 


—- linear regression - Eq. (2) 
+ model values 
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Fig. 12. Comparison of approximation accuracy relative to model values: 


H,= var, B = 21.42 m, CB = 0.56, GM = 0.8 m 
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— linear regression - Eq.(3) 
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Fig. 13. Comparison of approximation accuracy of the MSI relative to 
model values 
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Fig. 14. Comparison of approximation accuracy of vertical accelerations 
relative to model values 
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Fig. 15. Comparison of approximation accuracy of lateral accelerations 
relative to model values: at, ,,,,,,,- Model value of significant amplitudes 
of lateral accelerations, at „~ value of significant amplitudes of lateral 
accelerations, approximated by means of Eq. (5) 
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Fig. 16. Comparison of approximation accuracy of lateral accelerations 
relative to model values: GM = var, Lwl/B = 6.22, CM = 0.97, H= 3m 
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Fig. 17. Comparison of approximation accuracy of lateral accelerations 
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Fig. 18. Comparison of approximation accuracy of lateral accelerations 
relative to model values: XF = var, Lwl/B = 7.4, CM = 0.95, 
GM = 0.6m, H,=3m 


B [m] 
Fig. 19. Comparison of approximation accuracy of lateral accelerations 
relative to model values: B = var, Lwl/B = 7.0, CM = 0.97, 
GM — 0.8 m, H =3m 
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relative to model values: Lwl = var, Lwl/B = 6.6, CM = 0.95, GM = 0.6 m, 


Lwl [m] 


Fig. 20. Comparison of approximation accuracy of lateral accelerations 


H,=3m 


RANGE OF APPLICATIONS AND 
LIMITATIONS OF THE ELABORATED 
RELATIONS 


The elaborated design guidelines may find application to 


design calculations of the passenger-car ferries: 


of parameters and hull forms consistent with the ranges 


used for the model variants, 
intended to operate in the assumed service conditions. 


With taking into account the first group of the above 
mentioned limitations, the relations (1), (2), (3), (4) and (5) 
can be especially used for design calculations of passenger-car 


ferries of the following dimensions and parameters: 


3 


1 
2 
3 
4 
5. 
6. 
7 
8 
9 
1 


. LBd = 19 000 + 46 000 n 
. L/B =5.8+8.2 
. Bd =3+4.5 
CB =0.56+ 0.64 
CWL = 0.75 + 0.85 
LWL = 124 +258 m 
. B =19 =+ 33 m 
. CM =0.95 + 0.98 
. XF =55.2+117.6m 
0.Fw =2 100 -= 6 000 m’. 
CONCLUSIONS 


The elaborated design guidelines can be applied to 
preliminary ship design as target functions or technical 


limitations. 
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If recommended limit values of the presented sea-keeping 
qualities are known limit values of ship design parameters 
can be determined. For instance, to keep MSI value below 
20% is recommended according to [6]. And, by making use 
of Eq. (3) the limit area of waterplane can be calculated 
depending on the significant wave height Hs. 

The elaborated relations are characterized by different levels 
of complexity and approximation accuracy. Eqs. (1) and (5) 
are of the most complex structure (they were elaborated with 
the use of artificial neural networks), but simultaneously 
they show the highest approximation accuracy. 

All the relations are based on the data available in the 
preliminary stage of ship design process. 
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State equations in the mathematical model of 
dynamic behaviour of multihull floating unit 


Agnieszka Królicka, M. Sc. 
Gdansk University of Technology 


ABSTRACT 


This paper concerns dynamic behaviour of multihull floating unit of catamaran type exposed 
to excitations due to irregular sea waves. Dynamic analysis of multihull floating unit 
necessitates, in its initial stage, to determine physical model of the unit and next to assume 
an identified mathematical model. Correctly elaborated physical models should contain 
information on the basis of which a mathematical model could be built. Mathematical models 
describe mutual relations between crucial quantities which characterize a given system in 
time domain. The dynamic analysis of multihull unit was performed under assumption that 


the unit s model has been linear and exposed to action of irregular sea waves. Mathematical 
model of such dynamic system is represented by state equations. The formulated equations take into account 
encounter of head wave which generates symmetrical motions of the unit, i.e. surge, heave and pitch. For 
solving the equations the following three wave spectra were taken into consideration: 
- ISSC (International Ship Structures Congress) spectrum 
- Pierson-Moskowitz spectrum 
- Paszkiewicz spectrum. 


Keywords: dynamics of multihull floating unit, sea waves, stochastical processes, wave spectra 


INTRODUCTION 


Motion of a dynamic system can be generated by different 
external or internal factors. At mathematical modelling external 
excitation factors of the most significant effect on the system, 
are selected. Such external factors are usually called excitations. 
Response of the system to given excitations is mathematically 
characterized by a definite transformation called operator of 
a system. For a broad class of dynamic systems the relation 
between excitations and response is characterized by differential 
equations of motion. The equations can be linear or non-linear, 
of constant or varying coefficients, ordinary or differential, 
deterministic or stochastic ones. The mathematical models 
used for practical applications almost always necessitate to 
simplify the equations which form a given model. It amounts 
mainly to lowering their order by forming partial models and 
time — local ones with an assumed practical accuracy. Dynamic 
mechanical systems which represent floating objects are tightly 
associated with stochastic processes. State variables and 
input parameters of the models are of probabilistic character. 
Mathematical models of such systems are represented by 
sets of stochastic differential equations, and form sets of Itó 
Itequations. Multihull units such as catamarans and trimarans 
belong to complex, highly non-linear dynamic systems. If 
dynamic system model is of floating unit's linear system then 
the equations: 


6 
I; n; +B;n;+C;n = E(t) (1) 


i,j=l 


where: 
I=M+A- inertia matrix 
M — elements of matrix of generalized masses of 


structures 
A — elements of matrix of generalized hydrodynamic 
added masses 

B — hydrodynamic damping matrix 
C — hydrostatic stiffness matrix 
n — generalized displacement vector 
F(t) — vector of exciting forces and moments 
can be analyzed as a set of two mutually non-coupled groups 
of mutually coupled equations. It is assumed that the coupling 
takes place by linear and non-linear damping coefficients and 
hydrostatic stiffness coefficients. 

In these considerations the examined object is taken as 
a linear dynamic system of six degrees of freedom. 

Among them are the following: 
- translational oscillations: 

a) surging — n, 

b) swaying — n, 

c) heaving — n, 
- angular oscillations: 

d) rolling — 1, 

e) pitching — n, 

f) yawing — n, 


The catamaran-fixed coordinate frame is shown in Fig. 1. 
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Fig. 1. Schematic diagram of physical model of catamaran 


Local motions of the unit constitute its response to sea- 
wave-induced excitations. To the first group belong the 
equations which describe the symmetrical motions: n,, 1),; ns 
The second group is formed by the equations which describe 
the anti-symmetrical motions: n,, n,, Ne Constructional 
reasons (symmetry) of multihull floating units make it 
possible to analyze group-coupled motions of the object and 
in consequence to limit number of state variables which appear 
in the equations. 


SET OF STATE EQUATIONS 


State equations are one of the possible ways for representing 
mathematical model of a dynamic system. An alternative way 
to describe a dynamic system is transmittance in which the 
initial state is assumed equal to zero. 

The operator transmittance, otherwise called the transfer 
function G(s), is ratio of output signal Laplace transform and input 
signal Laplace transform of the system in zero initial conditions. 
The transmittance describes general features of a stationary 
linear system of one input and one output, independing on a kind 
of excitation. For the systems described by linear differential 
equations of constant coefficients the transmittance is rational 
function of the complex variable s, which can be represented by 
the quotient of two polynomials (2). 

To determine the set of state equations the following 
assumptions have been applied: 

1) resultant motion ofthe object on an irregular wave is formed 
by superposing its motions in regular waves 

2) only encounter of the head wave which generates only the 
motions: surging rj, heaving n, and pitching n,, is taken 
ito consideration 


3) it is assumed that the response of the object, in the form 
of the generalized wave-generated forces F, (i = 1, 3, 5), 
to excitation due to the wave of the height &(t), can be 
approximated by means of a system whose transmittance 
is of the form: 

F(s) _ bas? +b,s+b, 
ANI MOS dove f vm (2) 
E(s) = s’+a,s+a, 

ie 

FUE, PB) 

b bus b 03? bos), b (b,,. b, b,,), b 2(b,,, b 23? b, 


a (a, a 13? ays), a(a,,, a 23? 85). 


The relation (2) can be written by using the following set 
of state equations: 


f -F-h$ 
WM (3) 
f£, =F -hé - hé, 
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f=f thé 


; 4 
f, =-a, fi af +h » 
where: 
h, P: h, — constants defined by the coefficients ofthe equation 
T" p fi) £(£ A fas) 
h,(h 01? ho» hs). h (hh 13? his h,(h,,, h 23? h,;). 


To obtain random process of the irregular wave height 
&(t) the well-known energy (wave) spectra of the wave which 
encounters the object, can be used: 


e.g. Pierson - Moskowitz spectrum 
dee (0) = ) 
or ISSC spectrum 


The selected spectrum is approximated by means of the 

spectral density function in the form: 
g(o) = = (6 
o* —2vo,0° +03 

Then, the shape filter is introduced by making use of the 
following assumptions: 
1. In both the spectra, i.e.(5) and (6), their maxima appear at 

the same frequency and are of the same value 


2. fà, ()do = [g(0)do 
0 0 


3. The wave height processes are generated by the transfer 
function G(s) which constitutes the so-called shape filter 
if „white noise" is at input. 

The transfer function G(s) is ~ by the relation: 


Go) = — > 
B )- s? + asta, ) 
The relation (7) corresponds to the following set of state 
equations: 
g =È 
: (8) 
22 -6 -aW 
£ = 8 +aW 
n 9 
g, = 2,8 — & 8. — aa, W O) 
where: 


W - „white noise" 
&(t) - wave height process. 


The above given matrix (Tab. 1) can be written in the 
general form of Itó equation: 


X = AX + BW (10) 

From the matrix the following is obtained: 

Li =f, 

Li --ajf, —ay fy, 

f =f, 

fo =—a);f,; — auf, (11) 

fis = f, 

E m asf; E ass 

£j -g; +a, W 

£,—-—8,8,; - 4,8, - a,8,W 


Tab. 1. Matrix of state equations 


100000 be By Bo c, 0 0 
=|000000 =|b,, ba b, C;=|0 C4 C4 
000010 [ba b4 b, 0 C4 C55 
The parameters d, 4,, d, can be determined after finding Now o - value at which the extremum appears, is 
and solving the particular coefficients of damping, hydrostatic determined, i.e.: 
stiffness as well as inertia ones. Co’ 
g(0)= u^ TuLE (16) 
EXAMPLE SOLUTIONS FOR THE pred 
SELECTED WAVE SPECTRA le) [e 
as BI - B (17) 
1 ISSC wave spectrum From the condition that extremum is present in the same 
A B 12 point the following is obtained: 
d(0)= —; exp --7 (12) ic the first equation: 
© © ; 
where: = 1 
` 4B \4 = 4 
(=) =E* so Sr (18) 


he ZEJ 
Qi 


From the second assumption that: 


(13) 
T $0) =8(0) 
A= 173h? z 
i 
m 1 
The wave spectrum is approximated by means of the 0 = ^B hi 0 = E^ 
spectral density function of the following form: i i 
Co? 
8(9)- —5,— —3 3 (14) WERT 
Q —2vO,0 +0, the following is obtained: 
Let it be: i.e. the second eds 
2vo; =D E 


o, =E PE ENS E-D 
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and the third equation 

" 2 >. 
From the right-hand side of the Eq. (20): 

jes exp- Z Jo = afo” exp (Sjo (21) 


by substituting: 


0 


-40 do = dt > do = E 


the following is achieved: 


A [o^ expC Bi) + oat |= 
0 LS - 


(22) 

AT A - A 

=-— | exp(- Büdt - —e "| =-— 

4 J x ) 4B 0 4B 

The left-hand side of Eq. (20) is of the form: 
€ Co 

———— —— — do 23 
] ca*—ba* +a gi 


0 


in which ,,C” means ,,white noise", and the coefficients a, b, c 
are equal to: 


b - 2vo; 
acd Q4) 
c=l 


After solving the above given integral the final solution is 
achieved as follows: 


T Co? 
[—— do 
o —Do*~ +E 

-Cr 


20,4 VvV- e (Y? v 


By comparison ofthe left hand side and right hand side (the 
third equation) the following is obtained: 


i — - — = 
M 20, 6 4v? DE tWVtvv v 


This way the third, lacking equation for the unknowns C, 
v, ®,; has been achieved. 
The first assumption (Eq. (18)) yields the following: 


0 


(25) 


(26) 


4 4B 
E=a=0,=—B>0; + | (27) 
5 5 
The second equation (Eq. (19)) yields the relation: 
(28) 


SA(S Wi. C 
4B | 4B 2JE -D 
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By introducing the following substitution: 


D=b=2va; 
(29) 
E=0; 
the relation is achieved: 
4B \ 4B 203(1-v) 
The constants A and B are equal to: 
-4 
B=69 (= 
O, 
(31) 


-4 
A= vran( 2) 
O, 


From Eq. (30) is obtained the quantity C expressed as 
follows: 


2AU=% 
C= = 2 
o'e 
From Eq. (26) is determined the quantity C equal to: 


C= 250; ^ Jv? -1? + Gv? Di (33) 


47B 


Next the right hand sides of the above given equations are 
compared to each other to get the equation for determining 
the unknown ,,v". 


32) 


Vi cr T 
(34) 
1 4D +1 
Va um E 


Now the quantity L = 4.500406826 and its square 
L? = 20.2536616 is determined. 


Hence: 
v, =1.496512384 
Si (35) 
v, = 0002235192 
By substituting v, the following is obtained: 
a, =0.59607 (36) 
C = -0.068635565h; -@, (37) 
a, = VC z:0.261983903h, Jo, 
a, =0,42(l1-v)= 69 
=0.772010362a, -0.996506281= 
=0.7693131750, 
By substituting v, the following is obtained: 
a _ ou 0.054985159h;-@, (39) 
2502.977614 
a, = VC =0.234489144h, Jo, 
a, =0.7720103620, -0.891924669 = (40) 


= 0.6885750870, 


2. Pierson-Moskowitz wave spectrum 


A ( B 
1-4-4 (41) 
Te E (42) 
T, =0.92T, 
d (an) 1 ( a \ 
NT (E NETS E E 
5 
io =188.17 (43) 
Sr’ (0.927 | T, 


4 
1| 2% 
B=—|— | =496.1 
TG 
On the assumption that: L = 4.500406826 and 


v, = 1.496512384, v, = 0.602235192, value of the constant 
C is calculated: 


188.1787h20.710370, - 2.234507723 Hu 


C= 
= h? -œ -0.095826987 
a, = VC =0.309559344-h,- Jo, (45) 
0, = (S = 0.7103706810, (46) 
a, — 0, —0.5046265040? (47) 
a, = 6,4 2(1— v) = 0.9965062810,= m 
= 0.7078888440, 
3. Paszkiewicz wave spectrum 
A 2 0.812-10? 2m) _ =7.95 
T 
2 4 (49) 
B-0.23| C | =1126.8 
T, 
©, = j2- 0.7605261330, (50) 


a, =0; = 0.57840? (51) 


a,=®,y2(1-v) =0.996506280,= 
52 

=0.7578690670, ci 
C=hî-@,:1.9082547-10° (53) 

a, = VC -0.043683574-h,-Jo, | 69 


B ISSC spectrum 
JJ)" 
® Pierson-Moskowitz spectrum 


6 + m Paszkiewicz spectrum 
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Fig. 2. Distribution of values of the coefficients for the 
selected wave spectra 
Values of the functions for @, =1 and h, =1. 


CONCLUSIONS 


1. The coefficients a,, a,, a, which appear in the shape filter 
(7) were determined. They were determined with taking 
into consideration the three wave spectra: 

a) ISSC (International Ship Structures Congress) 
spectrum 

b) Pierson-Moscowitz spectrum 

c) Paszkiewicz spectrum 


2. The set of state equations which takes into account 
symmetrical motions of catamaran, was presented in 
the form of the matrix equation (given in Tab.2). They 
constitute the equations of dynamic behaviour of catamaran 
on irregular wave, called also Itó equations. The quantities 

, V, C obtained from Eqs. (18 + 20) made it possible to 
ae the coefficients a,, a,, a,. 

3. Further continuation of the subject matter in question is 
aimed at achieving a solution of the elaborated matrix 
equation. 


Tab. 2. Functional expressions for the coefficients of the selected wave spectra 


spectrum ot — -@,0.054985159 


Pierson-Moskowitz 
hg -œ -0.095826987 | hi 


Paszkiewicz 
-®, -1.9082547 - 10° 


0.234489144h, Jo, 


0.309559344-h, fo, 


0.043683574- h, -4 © 


0.6885750870, 


0.7078888440, 
0.5046265040; 


0.7578690670, 
0.57840; 
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Experimental verification of two-parametic 
models of fatigue characteristics by using 
the tests of SS5J0 steel as an example 


Bogdan Ligaj, Ph. D. 
Grzegorz Szala, Ph. D. 
University of Technology and Life Science, Bydgoszcz 


ABSTRACT 


This paper presents experimental verification of models oftwo-parametric fatigue characteristics N(0.,, a...) 

on the basis of results of fatigue tests of specimens made of S355J0 steel. On the basis of comparative 

analysis of results of computations and fatigue life tests of the specimens, performed under constant amplitude 

sinusoidal loads of different cycle asymmetry ratios, the most versatile models were distinguished. The 
tests in question were carried out in the high-cycle fatigue range. 


Keywords: two-parametic models, fatigue characteristics, tests of S55JO steel, 
computations, fatigue life tests, high-cycle fatigue range 


INTRODUCTION 


Random loads commonly occur in service conditions 
of machines and devices including sea-going ships, ocean 
engineering objects, land road vehicles, aircrafts, heavy 
machines etc. From the point of view of mathematical theory 
of stochastic processes operational loads are non-stationary 
stochastic processes [1, 2], that makes direct application of 
the theory in question to calculating fatigue life of structural 
elements, difficult. 

In practice fatigue life of structural elements is calculated 
in accordance with the schematic diagram shown in Fig. 1. 
As results from it, the following should be known to carry 
out fatigue calculations: service load spectrum and fatigue 
characteristics, as well as to assume an appropriate hypothesis 
of cumulating fatigue damage [3]. 


SERVICE LOAD SPECTRUM 


HIPOTHESIS OF CUMULATING FATIGUE 


MATERIAL CYCLING PROPERTIES 


(material fatigue characterstics) 


DAMAGES 


CALCULATED FATIGUE LIFE OF 
STRUCTURAL ELEMENT 


Fig. I. Schematic diagram of calculation process of structural element 
fatigue life 


In fatigue life calculations random run of service loads is 
usually substituted by a set of sinusoidal cycles determined 
in accordance with appropriate methods [4]. Their spectra 
can be described by means of distribution of load amplitude 
substitute values as in the case of broad-band random loads 


-by distribution of S, - amplitude values and S, - mean values 
of cycles or S... - minimum values and S|. - maximum 
values. Fig. 2 shows a scheme of possible cases of sinusoidal 
cycles having different parameters. The sinusoidal cycle is 
unambiguously described by the following parameters: S 
Sat = VEorS, S. f 

Influence ofthe load frequency fon fatigue life of structural 
elements is usually low and neglected in calculations of many 
elements. Location of a load cycle within the reference system 
(S, o 9, OF (S, S,) is determined by the cycle asymmetry 
ratio R = S,./S.., or sometimes — the load steadiness ratio 
X7 S JS, 

Load spectra are then elaborated in one ofthe two reference 
systems: (S... S...) [5] or (Sp S.) [4, 6]. In such cases 
appropriate fatigue characteristics is that characterized by two 
parameters: N(S,, S.) or N(S,... Spad) 

In this paper the problem is discussed of experimental 
verification of selected two-parametric fatigue characteristics, 
performed on the basis of test results of specimens made of 
8355J0 steel. 


min? 


DESCRIPTION OF SELECTED 
TWO-PARAMETRIC FATIGUE 
CHARACTERISTICS 


Concept of description of limit stress diagrams (fatigue limit 
range) with taking into account mean value of cycle amplitude 
was presented in [7]. The descriptions of range of limited 
fatigue life were presented in [8, 9, 10, 11, 12, 13]. 

In Tab. 1 are collected the formulae describing the above 
mentioned characteristics together with references to their 
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Oscillating load, 


n 
S max 


Pulsating 
compression load 


Static compresion 
load, R = 1 


Pulsating tension 
load, R = 0 


Sm 


Static tension 
load, R = I 


Variability of R within the 


intervals: 
1) O<R<l| 
2) -1<R 
3) -0<R<-l AS 
4) 1<R<+0 


Fig. 2. Plane of variability of sinusoidal load parameters 


literature sources. The formulae (4), (5) and (10) concern 
description of limit stresses within the range of fatigue limit 
called also the unlimited fatigue life (FL) (N—00) for steel. In 
practice it is assumed that it concerns the cycle number N 7 N,, 
and N, is the cycle number corresponding with the slope change 
point on Wóhler fatigue diagram. The remaining formulae: (1), 
(2), (6), (7), (8), (9), (11) and (12) deal with description within 
the high-cycle fatigue range (HCF), called also the limited 
fatigue life. The range is described by the conditions: 


Sa Sa TS SR - for elastic-brittle materials 


S =S, Sa ER -for elastic-plastic materials. 


max 


In the first three positions of Tab 1. are given the relations 
corresponding with limit stress diagrams in the reference system 
(S S,) for unlimited fatigue life (N00). The simplest form of 
Goodman's relation, i.e. the formula (1), corresponds with straight 
line in the above mentioned reference system. Gerber’s limit 
stress diagram, i.e. the formula (2), is of parabolic form. Haigh’s 
diagram, i.e. the formula (3), is based on elliptic relationship. 

The two-parametric fatigue characteristics analyzed in this 
paper constitute a generalized form of the above mentioned 
relations used for description of limit stress diagrams in the 
range of high-cycle fatigue (HCF), i.e. the limited fatigue life 
range. The relations are here shortly called models. The oldest is 
Heywood’s model (H) described in the subject-matter literature 
[4, 8], which is used, in this paper, for comparision with the 
remaining models. The model I, an expanded form of Haigh’s 
relation, consists of two straight line sections; the first of them 
corresponds with the variability range of the cycle asymmetry 
ratio R between 0 and 1, the second — between 0 and -oc. 


The diagram in the broken-line form crosses the three points: 
(S, — 0; R or R), (7 RS ^ RY) and (S „= 0; RY). The model 
II, a generalized form of Goodman ’s relation, covers, by single 
relation, the variability range of R from 1 to -c0. The diagram 
lines of the model cross the two points: (S, = 0; R_ or R,) and 
=0; RT). 
The "model III is an expansion of Gerber's relation [the 
formula (2)], and the model IV is an expansion of Haigh's 
relation [the formula (3)]. In both the above mentioned models 
the diagram lines cross the points: (S, = 0; R, or R,) and 
(S, 2 0; RN). 

The model V which connects features of the models: H, I 
and III, is based on parabolic relation where the parable crosses 
the three characteristic points as in the model I. 

The last model, acc. Bototin, is a transformation of the 
model I into the reference system (S, ., S...) and hence it has 


min? ~ max 


been not subjected to any detail analysis. 


EXAMPLE TWO-PARAMETRIC FATIGUE 
CHARACTERISTICS FOR S355J0 STEEL 


m 


Strength properties of S355J0 steel under monotonic loads 
is given in Tab. 2, and its cyclic properties — in Tab. 3. 


Tab. 2. Static strength properties of $355J0 steel 
Static properties of S355J0 steel 


Standard deviation] 84 | 71 | 1306 [099] 09 


Tab. 3. Cyclic mechanical strength properties of $355J0 steel 


Pulsating load "A 
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Tab. 1. Collection of the formulae which describe fatigue characteristics 


Name and form of formula 


Fatigue diagram acc. Goodman 


Fatigue diagram acc. Heywood 
S S 
ps ( -$n Ja, +1(1-A0)] 


where: Ao andy - constants determined with the use of smooth 
specimens: 
- for steel: 


S | se) 1+0.0038(logN)* 
Y 248.|, Ay =———— 
1+0.008(logN) 


- for AI alloys: 
" 0.003 1(log N)* 
|... 1+0.003(logN)* 
7 F strane | '^* 1+0.003I(logN)* 
225 
Fatigue diagram acc. Szala (the model I) 


N,R™ 
(s, * WwSm 7" 


NaN RaRa 8, 78,2 | " 
^| SR, w,) 


where: 
de cdi (4-4) 
Wy 22€, C ™NS™ ™ 1 (7) 
Formulae (6) and (7) presented in another form: 


1 
-y Sm, Raf No \™ 
YR RUN 


m 


[ 
[12] 


Application range 
Variability range R 


Literat. 
source 


t e ges I I Nl NND 


Range of fatigue limit 
(FL) 
Sa <R, 


Range of fatigue limit FL 
FL i 
S SR, 
Range of fatigue limit FL 
FL 
-æ <R <1.0 [14] 


7] 
[4] 
[8] 

9] 


[14] 


Range of HCF for: 


- brittle materials 
Six =S, +S, <R, 


- elastic-plastic materials 
Six =S, +S, SR, 
—o«R «1.0 


Range of HCF: 


-o<R<0 


0<R<1.0 


-0<R<0 


0<R<1.0 


POLISH MARITIME RESEARCH, No 1/2010 


41 


Tab. 1. Collection of the formulae which describe fatigue characteristics 


No. Name and form of formula Application range Literat. 
Variability range R source 


Range of HCF: 


vt Dm -æ <R <1.0 


-œ < R <1.0 


Range of HCF: 


-œ <R <1.0 


Range of HCF: 


-æ <R <1.0 


Fatigue diagram acc. Bototin Range of HCF: 
NS NR? — F id nP 
[A+ Y)S max - (17 W)S min 1 "AC SK 
_ 2" N R 
[O — y')S max — (A+ WS min ]"? 
where: y and y'- material constants 
factors of material sensitivity to cycle asymme 
Comments and basic formulae: 
Wöhler diagram for R = -1: 
- general form: S? -N =C, 


for 


- logarithmic form: logS, = ioe Ne logC, 
m, 


Wöhler diagram for R = 0: 
- general form: S™ -N =C 


- logarithmic form: logS,.. = d gN +logC 
m 


Saa Sat 
: R = min — m a 
cycle asymmetry ratio S. (S. 48.) n S 
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Knowledge of mechanical properties of S355J0 steel makes 
it possible to determine two-parametric fatigue characteristics 
complying with the models presented in Tab. 1. The models 
specified in items: 1, 5, 6, 7, 9 and 10 have been taken into 
consideration for further analyses. 


Heywood’s model 


The Heywood’s model (item. 4, Tab. 1) was selected 
to exemplify a way of determining two-parametric fatigue 
characteristics, their graphical form and characteristic features 
of diagrams. Such diagrams are elaborated in the reference 
system (S. , S), in the form of contour line diagrams. Particular 
contour lines are of constant fatigue life N. Hence when an 
appropriate value N is assumed a corresponding contour line 
can be determined on the basis ofthe formulae (4) by calculating 
S, - value for varying values of R. For subsequent values of N 
subsequent contour lines are determined the same way. 

The contour line diagram for S355J0 steel, according to 
Heywood's formulae, is shown in Fig. 3 for N = 10°, 104, 105, 
10° and 10°. 


Smax = const 


Smin = const 


7 
Pa S, = const f 


4 
4 


-700 -600 -500 -400 -300 -200 -100 


= Amplitude Sa, [MPa] 


As results from Fig. 3 varying values of S, - stresses take 
place within the range (-R, ; +R _), and those of S, - amplitudes 
within the range (0; R_,). To the diagram are also introduced 
example lines corresponding with the directions of constant 
values of: 

- cycle asymmetry ratio (R = const) 
- cycle maximum stresses (S... = const) 
- cycle minimum stresses (S... = const) 
- amplitude of stresses (S, = const) 
- cycle mean stresses (S„ — const). 


Similar diagrams were elaborated on the basis of calculation 
results according to the models: I - (see Fig. 4), II - (Fig. 5), III 
- (Fig. 6), IV - (Fig. 7) and V - (Fig. 8). 

From a general assessment of runs of diagram's constant 
-value lines for the assumed values of N, significant differences 
in the runs for the particular models can be observed. 

For practical applications of the discussed models 
of two-parametric fatigue characteristics conformity of 
results of calculations and verifying tests is of fundamental 
importance. 


I 
I 
I 
I 
I 
I 

n 
I 
M. Sm = const 

I 


100 200 300 400 500 600 700 


Mean value Sm, [MPa] 
Fig. 3. Two-parametric fatigue diagram acc. Heywood model (formula 4) for S355J0 steel 


Run — ——7100 
—»—> > 


Amplitude,Sa, [MPa] 


-700 


— o4 
-600 -500 -400 -300 -200 -100 0 


100 200 300 400 
Mean value Sm, [MPa] 


Fig. 4. Two-parametric fatigue diagram acc. Szala’s model (formula 5 and 6) for $355J0 steel 


500 600 700 
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-700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 
Mean value Sm, [MPa] 


Fig. 5. Two-parametric fatigue diagram acc. Szala s model (formula 8) for $355J0 steel 


Amplitude Sa, [MPa] 


-700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 
Mean value Sm, [MPa] 


Fig. 6. Two-parametric fatigue diagram acc. Lipski s model (formula 9) for $355J0 steel 


Rm ——700 — Amplitude Sa, [MPa] 
— 


= -1.25 —| 


-700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 
Mean value Sm, [MPa] 


Fig. 7. Two-parametric fatigue diagram acc. Lipski s model IV (formula 10) for $355J0 steel 
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-700 


-600 -500 -400 -300 -200 -100 


Rm ——700 A Amplitude Sas [MPa] 
———————Ó- 


100 200 300 400 500 600 700 


Mean value Sm, [MPa] 
Fig. 8. Two-parametric fatigue diagram acc. Lipskis model V(formula 11) for S355J0 steel 


VERIFYING TESTS 


Quantitative analysis of conformity of results of calculations 
and verifying tests can be performed on the basis of comparison 
of Wóhler fatigue diagrams for definite values of the cycle 
asymmetry ratio R. In this research the following values of R 
were assumed: R = 0,5; 0; -0,5; -1,0; -1,25; -2,0 and -3,0. The 
values cover to a large extent the parameter variability area of 
sinusoidal cycles, shown in Fig. 2. 

Results of the experimental tests and Wóhler fatigue 
diagrams in the reference system S, (N), determined on their 
basis for the above mentioned values of the cycle asymmetry 
ratio R, are given in Fig. 9. The diagrams in the bi-logarithmic 
form are described by the equations (16): 


logS, = alogN + b 


Values of the parameter a and b are given in Tab. 4. 

Data contained in Tab. 4 make it possible to present 
results of the tests in the form of a contour - line diagram 
of two-parameter fatigue characteristics (Fig. 10) similar to 
the diagrams based on results of calculations carried out in 
compliance with the analysed models (Fig. 3 through 8). 


(16) 


Stress S, [MPa] 


200 


10° 10° 10° 10° 10° 107 
Number of cycles, [N] 


Fig. 9. Results of the fatigue tests and Wöhler fatigue diagrams determined 
on their basis, for the selected values of the cycle asymmetry ratio R 


Amplitude Sa, [MPa] 


-700 


! oe 
-600 -500 -400 -300 -200 -100 0 


100 200 300 400 500 600 700 


Mean value Sm, [MPa] 


Fig. 10. Two-parametric fatigue diagram elaborated on the basis of experimental data for $355J0 steel 
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Tab. 4. Values of the directional coefficient a and the free term b in the 
equations describing Wöhler fatigue diagram, for the assumed values 
of the cycle asymmetry ratio R 


Cycle asymmetry ratio R 


-0.5 -1.0 -1.25 


Number of cycles [N] 


Fig. 11. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R = -3 


Stress S. [MPa] 
a 


Number of cycles [N] 


Fig. 12. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R = -2 
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In the subsequent figures, Fig. 11 through 17, comparison is 
presented of calculation results obtained by using the selected 
models of two-parametric fatigue characteristics for the particular 
values of the cycleasymmetryratio R, withexperimental testresults. 
In Fig. 12 through 16 the calculation results are presented on 
the background of experimental diagrams. The shadowed zones 
indicate breadth of experimental data scatter bands. Moreover, 
in the diagrams a maximum value of S... complying with the 
condition Sax < R is given. 


X 


Stress S, [MPa] 


VENERE UN 


300 


200 


- 10° 10° 10° 10° 10° 10” 
Number of cycles [N] 


Fig. 13. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R = -1.25 


Stress S, [MPa] 


Sa max = 678 MPa 


300 


200 


150 


10° 10° 10° 10° 10° 10’ 
Number of cycles [N] 


Fig. 14. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R = -1.0 


In Fig. 11 (for R = -3) and in Fig. 17 (for R = 0.5) the 
experimental fatigue diagrams were determined by extrapolating 
the test results for the range of R (-2,0; 0). 

From overall assessment of the data shown in Fig. 1 through 
17 result differences between values of the calculated stress 
amplitudes S, and those obtained from tests, S , 

depending on an assumed calculation model, assumed 
value of fatigue life N and values of the cycle asymmetry 
ratio R. 


Stress S, [MPa] 


E BEE NE 
Sa max = 508.5MPa 


—#- Haywood's model 
==" Experiment 


Number of cycles [N] 


Fig. 15. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R = -0.5 


Stress S, [MPa] 


© Modell 

A Model II 

© Model III 

9 Model IV 

* Model V 

E- Haywood's model 
== Experiment 


Number of cycles [N] 


Fig. 16. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R — 0 


Stress S, [MPa] 
"T Ls 
© Modell 
A Modelli 
400 B Model III 
© Model IV 


* Model V 


E- Haywood's model 


300 


Sa max = 169.5MPa 


-— FÉ 

—— "JÉ EE — 
10° 10° 10° 10° 10° 10° 
Number of cycles [N] 


Fig. 17. Results of fatigue calculations according to the assumed models 
of two-parametric fatigue characteristics for R =0.5 


ANALYSIS OF THE CALCULATION 
AND TEST RESULTS 


Differences between results of calculations of S, and results 
of tests of S... related to amplitude values depend on a model 
of assumed values of the fatigue life N and cycle asymmetry 
ratio R, taken for calculations. 

The relative differences, 6, between results of calculations 
and tests, determined according to the formula (17): 


B 
Bes 


—_ —..10096 (17) 
are presented in the form of bar diagrams in the successive 
figures as follows: for the Heywood's model — in Fig.18, 
model I - Fig.19, model II - Fig. 20, model III - Fig. 21, model 
IV- Fig. 22 and model V - Fig. 23. Within the entire variability 
range of the fatigue life N and cycle asymmetry ratio R the 
largest differences appear in the Heywood's model (H), and in 
the extreme cases they reach values from +49.32% to -23.21%. 
The largest differences occur for high values of fatigue life 
(N = 10° + 10’). The interval corresponds with stress values 
observed in the cycles ofa significant share in service load spectra 
of structural elements. The fact greatly affects fatigue calculations 
whose aim is to assess serviceability of structural elements. 

The above mentioned differences result from that the 
Heywood’s model was based on the data of material strength 
properties obtained from monotonic tension tests. 

In the remaining models material properties are represented 
by the fatigue life determined either from the test under 
oscillating load (R = -1) — in the case of the models: II, III, IV 
and V, or the test under oscillating load (R = -1) and pulsating 
load (R = 0) — in the case of the model I. Therefore the 
calculation results according to the above mentioned models 
for values of R = -1, and additionally to the model I for R = 0, 
are in full conformity with the test results, that is obvious and 
in particular diagrams marked: 6 = 0. 
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Out of all the analyzed models for the entire range of fatigue 
life N and cycle asymmetry ratio R the calculations according 
to the models: I, III and IV showed the greatest conformity 


with the experiment. 


48 


b 
> 


t2 
© 


Difference in amplitude values [%] 
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10° 10° 10° 10° 10° 107 
Number of cycles [N] 


Fig. 18. Diagram showing difference in amplitude values determined 
with the use of Heywood s model 
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Fig. 19. Diagram showing difference in amplitude values determined 
for the model I 
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Fig. 20. Diagram showing difference in amplitude values determined 
Jor the model II 
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In the case of the model I the calculations for R = 0.5 


constitutes an exception, whose discrepancies from test results 
reached about 65% for N = 10°. 
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Fig. 21. Diagram showing difference in amplitude values determined 
for the model III 
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Fig. 22. Diagram showing difference in amplitude values determined 
for the model IV 
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Fig. 23. Diagram showing difference in amplitude values determined 
for the model V 


The model I is sensitive to magnitude of the factor of 
material sensitivity to cycle asymmetry, y,. In the case of the 
results shown in Fig. 19 the calculations were carried out with 
the use of y, - value determined from Eq. (7), Tab. 1. Data 
for the calculations were taken from Wóhler diagrams for the 
load variability ratios: R = -1 and R = 0. The large scatter of 
fatigue test results on the basis of which the diagrams were 
determined, makes accuracy of determination of the factor yy 
from Eq. (7), low. 

By modifying the model I with regard to magnitude of 
the factor y,, conformity between calculations performed for 
the model and appropriate experimental data can be greatly 
improved. 


CONCLUSIONS 


From the above presented analysis the following results: 


a) conformity between calculation results and fatigue test 
results greatly depends on assumed two-parametric fatigue 
characteristics 


b) within the entire range of the analysis the highest conformity 
was obtained for the models: I, III, IV, and the lowest - for 
Heywood’s model 


c) the low conformity of calculations according to Heywood’s 
model results from the assumption of the tensile strength R 
as a quantity characterizing properties of material subjected 
to cyclic loads 


d) the higher conformity of calculation results according to the 
remaining models: I through V, results from the assumption 
of Wohler fatigue diagram parameters as characteristics of 
material mechanical properties; such data for a variety of 
materials can be obtained from tests or literature sources 
of different kind, e.g. [16, 17] 


e) for preliminary calculations, e.g. structural design process 
in which precise data on material cyclic properties and 
service loads are not yet known, the simple models I and 
II are recommended 


f) for further considerations an analysis of influence of 
the differences in two-parametric fatigue characteristics 
described in this paper, on calculation results of fatigue 
life in conditions of service load spectra, is of a great 
importance 


g) as results from literature data, especially those contained 
in [7], the above described models can yield different 
conformity with results of tests of other materials, the 
opinion is confirmed by the paper [12] devoted to tests of 
aircraft Al-alloys 


h) another important problem is determining range of 
application of the described models to high-cycle fatigue 
(HCF), which requires determining an area of two- 
parametric fatigue characteristics in compliance with 
the assumed criterion [17]; as the criterion S|. < R 


x e? 


commonly met in the subject-matter literature, is considered 
approximate and ineffective in some cases 


1) a seperate problem important from the point of view of 
practical applications is description of two-parametric 
fatigue characteristics for notched elements. 


The problems specified in the points: f), g), h), 1) will be 
discussed in next papers. 


NOMENCLATURE 

A — elongation [% ] 

C — constant in the formula describing Wóhler 
fatigue diagram for off-zero pulsating load 
(R= 0) 

C, — constant in the formula describing Wóhler 


fatigue diagram for oscillating load (R = -1) 

— cycle number - general notation (fatigue 
life) 

— cycle number — fatigue life corresponding 
with fatigue limit 

— cycle asymmetry ratio 

— material yield point [MPa] 

— fatigue limit - general notation [MPa] 

— material tensile strength [MPa] 

— fatigue limit under pulsating load (R = 0) for 
N, cycle number, [MPa] 

— fatigue strength under sinusoidal pulsating 
load (R = 0) for N cycle number, [MPa] 

R — fatigue limit under oscillating load (R = -1) 

for N, cycle number, [MPa] 
RN — fatigue strength under sinusoidal oscillating 
load (R = -1) for N cycle number, [MPa] 

— specimen stress — general notation, [MPa] 
max” Sin) © Sinusoidal cycle stress amplitude [MPa ] 
—0,5(S. + Snin) - mean sinusoidal cycle stress [MPa ] 

" — maximum sinusoidal cycle stress [MPa ] 

— minimum sinusoidal cycle stress [MPa ] 

— contraction [%] 

— exponent in formula describing Wohler 

fatigue diagram for pulsating load (R = 0) 

m, — exponent in formula describing Wohler 
fatigue diagram for oscillating load (R = -1) 

y — factor of material sensitivity to cycle 
asymmetry, for N=N, 

Wy — factor of material sensitivity to cycle 
asymmetry, for N # N,. 


2 2 
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min max 
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Comparing N3-60 cascade exit angles obtained 
from tunnel measurements and numerical 
simulations 


Biernacki Rafat, Ph. D. 

Gdansk University of Technology 
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ABSTRACT 


The article compares the results of the measurements performed in the aerodynamic tunnel owned by the 
Czestochowa University of Technology and the numerical simulation done using the code FLUENT for 
a plane cascade of N3-60 profiles. The comparison study aims at assessing differences between the exit 
angles measured experimentally and those obtained from numerical calculations. Variable parameters in the 
measurements and the calculations were the relative pitch and the profile stagger angle. The measurements 
were performed in the flow of air. The same medium and having the same inlet parameters was assumed in 
the calculations. The calculations were performed in two dimensions for the compressible flow. The k-e RNG 
turbulence model was used to complete the equation system. Comparing the measured and calculated results 
provides opportunities for assessing the range of differences between the experiment and the simulation. 
The obtained results and formulated conclusions will find the application not only in the power industry, 
but also in the shipbuilding industry for analysing the operation of steam and gas turbines used as main 
and auxiliary propulsion systems on both merchant vessels and battle ships. 


Keywords: numerical calculations, profile characteristics, gas turbines, steam turbines 


GEOMETRY OF N3 PROFILE and the relative pitch t/b we can calculate the effective cascade 


exit angle, see the diagram in Fig. 3. 


The N3-60 profile is a basic profile used in impulse turbines 
produced in Poland. In the eighties, a series of tests of those 
profiles were performed in the aerodynamic tunnel owned by 
the Czestochowa University of Technology. Variable parameters 
in these tests were the profile stagger angle o, and the relative 
pitch t/b. The tests aimed at preparing basic flow characteristics 
of the profiles and cascades to use them in the turbine design 
process. An isolated profile N3 is shown in Fig. 1. 


Fig. 1. N3 profile contour lines [1] 


Basic geometrical parameters which define the N3-60 


profile are given in Fig. 2. For the changing stagger angle a, Fig. 2.:Geometrical parameters of N3-60 profile. Ref 
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0.5 0.6 0.7 0.8 


0.9 1.0 1.1 1.2 t/b 


Fig. 3. Effective N3-60 cascade exit angle [1] 


DISCUSSION OF THE EXISTING RESULTS 
OF N3 PROFILE MEASUREMENTS 


The existing results of N3 profile measurements are 
published in [1]. The ranges of parameters changing in the 
profile tests were the following: 

* relative pitch t/b - 0.5 to 1.2 

inlet flow angle - 45°, 75°, 90°, 120° 
cascade stagger angle - 34.6° to 48.6° 
inlet flow velocity - Ma=0.2 

resultant Reynolds number - 0.38x106 


The measured results, having the form of velocity 
distributions, were integrated and averaged to obtain one- 
dimensional characteristics of the cascade for the assumed 
profile stagger angle and pitch. The quantities determined in 
those measurements included the profile loss, the trailing edge 
loss, and the exit flow angle a,. Fig. 3 shows the flow exit angle 
as a function of the relative pitch for the inlet flow angle a, equal 
to 90°. The basic goal of this study was to verify the differences 
between the real results obtained from the measurements and 
those calculated using the numerical model of the cascade. 


AUTOMATION OF THE GEOMETRY 
GENERATION PROCESS 


The profile co-ordinates were entered to the Excel 
calculation sheet. The entered data included fifty pairs of points 
in the Cartesian coordinate system which defined the upper and 
lower profile curve. The next step was to calculate the camber 
line, defined as: 

x = X, =X d 


y, 120 t Ya) 


For an arbitrary N3 profile stagger angle, the algorithm 
rotated the profile and its camber line using the following 
scheme: 
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x’ = xcos(a) — ysin(a) 
y’ = xsin(a) + ycos(a) 


After the rotation X, # X, £ X. The passage axis was 
generated by shifting parallelly the camber line up and down 
by half of the assumed cascade pitch, and then extrapolated 
beyond the profile as the straight lines extending from 
x=-25 tox = 100. The input data for the geometry generation 
algorithm comprised a list of points which defined the profile 
in standard position, complemented by the rotation angle and 
the cascade pitch. Based on these data, the algorithm generated 
a script which introduced to the code Rhino a sequence of 
points describing the passage geometry after rotation and 
transformation. This method was used for creating a series of 
passage geometries taking into account different profile stagger 
angles and cascade pitches. Fig. 4 shows a sample geometry 
before and after the transformation for the stagger angle 
a = 44.58? and the relative pitch t/b = 0.6. 


after rotation 
——— Passage axis 


Fig. 4. Profile transformation and creation of the passage axis 


The oriented N3 profile, mapped using the code Rhinoceros, 
was used as a basis for generating the calculation grid. The 
2D calculation domain was bordered in front by a vertical 
inlet section, and at rear by a vertical exit section. From the 
top and the bottom the calculation domain was bordered by 
the axis of the blade-to-blade passage, treated as the border 
of periodical type. 


Profile N3-60 
alpha 44.58° 
vb = 0.6 


[s 


Fig. 5. Geometry prepared for grid generation 


ASSUMPTIONS ADOPTED WHEN 
GENERATING THE CALCULATION GRID 


Inall cases, the geometry presented in the previous Section 
has made the basis for generating the calculation grid. The 
upper and lower boundaries of the calculation domain (see 
Fig. 5), which define the flow passage axis, were considered 
the periodic boundaries. The calculation domain was divided 
into two subdomains — over and below the plate, which made 
it possible to generate the boundary layer on the surface of 
the profile washed round by the flow. The calculation grids 
of structural type and dimensions 160 x 450 elements were 
generated automatically by the code, which made it possible 
to generate an optimised structure for each examined geometry 
variant. The generated grid was optimised with respect to 
the orthogonality of the gridlines. One of basic criteria for 
assessing the quality ofthe created grid was the minimum angle 
between the lines which limit the grid element. Tab.1 collects 
the values of the minimal angles for different calculation 
variants. 


Tab. 1. Minimal calculation grid angles 


Profile stagger angle a, 
40.58 44.58 
2.70? 3.55° 
3.44° 5.49° 
3.15? 6.89? 


36.58 
3.06? 
3.64? 
3.91? 


48.58 
4.47? 
6.44? 
7.78? 


Fig. 6 presents a sample calculation grid generated for the 
variant characterised by parameters: a, = 44.58°, t/b = 0.6. 


Fig. 6. Sample calculation grid for N3-60 profile 
RANGE OF CALCULATIONS 


The calculation grid similar to that shown in the precious 
Section was generated for each individual pair of profile stagger 
angle a, and pitch t/b. The calculations done using these grids 
made use of parameters, the range of which is given in Tab. 2: 


Tab. 2. Range of parameters used in calculations 


Profile stagger angle o, 
40.58 | 44.58 


36.58 48.58 


The calculations were performed using the code Fluent [2]. 
The boundary conditions for all variants were taken from the 
experimental data. The assumed quantities included: 

- mass flow rate through a single passage 

- pressure at flow area exit 

- type of working medium (the experiment was carried out 
for air) 

- inlet temperature of the working medium 

- model of fluid (viscous, compressible) 

- turbulence model (k-e RNG) 

- values of turbulence parameters k and e at passage inlet. 


Part ofthe above parameters was constant for all calculated 
variants. This list includes: total inlet pressure assumed as 
equal to 100 kPa, inlet temperature of the working medium 
-300 K, and the model and parameters of the working medium. 
Boundary conditions which changed from variant to variant 
are given in Tab. 3. 


Tab. 3. Variable boundary conditions 
m [kg/s] 

2.480 

3.307 


k [m?/s?] 
21.627 
21.627 
21.627 


£ [m?/s?] 
2112.7 
1584.5 
1267.6 


Cascade 0). 


pitch 0.8 
t/b d 4.133 


The values assumed in the calculations imitated the flow 
conditions recorded in the experiment. 
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SAMPLE RESULTS OF NUMERICAL == 
CALCULATIONS 2805001 


The obtained results of the calculations have made it at I 
possible to assess the effect of the distance of the test section 2200401 
from the trailing edge on the values ofthe exit angle a,. As could 
be expected, the calculated angle did not change much after 1508-01 
leaving the profile cascade. Changes of this angle as a function 1.50601 
of the distance from the blade trailing edge are shown in Fig. bed 
7, for the selected sample pitch t/b — 0.6. When comparing jm ON P 
with the results of the measurements, it was assumed that the 8 00e«00 ) 


distance from the trailing edge was equal to 5 mm. 5 00e:00 
254 


c) 


Fig. 10. Results of calculations for N3-60 profile, a, = 44.58°, t/b = 0.6; 
a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 
d) changes of exit angle a, along datum +5 mm 


3.00e+01 = \ | = 
E HER \ 
= ie — 9 
0 I \ 
pa \ 
2.10e*01 \ 
— alfay 36.58 206501 b) 
5 — alfay 40.58 1800 
— alfay 44.58 1900-01 a) 
— alfay 48.58 Ex 
04 i 1.20e+01 
0 5 10 15 1100001 
Distance behind trailing edge [mm] g00e:08 d) 
Fig. 7. The angle a, vs. the distance from the trailing edge & 008100 
5.00e+00 
s . 4.00e+00 
The next figures show the streamline and exit angle same 
distributions in the calculation domain, as well as the exit angle 0000+00 | 
distributions along the traverse line x = +5 mm downstream Fig. 11. Results of calculations for profile N3-60, a, = 48.58°, t/b = 0.6; 
of the trailing edge. a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 
~ x d) changes of exit angle a, along datum +5 mm 
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Fig. 8. Results of calculations for profile N3-60, a, = 36.58°, t/b = 0.6; : " : - sa - 
a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); Fig. 12. Results of calculations fc or profile N3-60, a, = 36.58", t/b = 0.8; 
d) changes of exit angle a, along datum +5 mm a) streamlines; b) streamlines m enlarged; c) exit angles (scale on the left); 
— d) changes of exit angle a, along datum +5 mm 
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Fig. 9. Results of calculations for profile N3-60, a , = 40.58°, t/b = 0.6; Fig. 13. Results of calculations for profile N3-60, a, = 40.58°, t/b = 0.8; 


a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 


a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 
d) changes of exit angle a, along datum +5 mm 


d) changes of exit angle a, along datum +5 mm 
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Fig. 14. Results of calculations for profile N3-60, a, = 44.58°, t/b = 0.8; Fig. 18. Results of calculations for profile N3-60, a = 44.58°, t/b = 1.0; 
a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 
d) changes of exit angle a, along datum +5 mm d) changes of exit angle a, along datum +5 mm 
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Fig. 15. Results of calculations for profile N3-60, a, = 48.58°, t/b = 0.8; . : 
a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); Fig. 19. Results of calculations for profile N3-60, a, = 48.58°, t/b = 1.0; 
d) changes of exit angle a, along datum +5 mm a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the left); 
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Fig. 16. Results of calculations for profile N3-60, a, = 36.58°, t/b = 1.0; 251- 
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Fig. 21. Comparing exit angles a, obtained from measurement 
and calculations as a function of pitch t/b and stagger angle a, 


Fig. 17. Results of calculations for profile N3-60, a,= 40.58°, t/b = 1.0; From the point of view of impulse turbine flow system 


a) streamlines; b) streamlines — enlarged; c) exit angles (scale on the lef); calculations these difference are of high significance. For the 
d) changes of exit angle a, along datum +5 mm 
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i œ [0] 


smallest stagger angles (about 36°) this difference approximately 
equals to 3°. Taking into account that the exit angles in this area 
are within 12°-15°, the above difference corresponds to an error 
of an order of 20-30 %, which then leads to a similar error in 
assessing the mass flow rate though the turbine flow system. 
To better illustrate this problem, Fig. 22 shows the differences 
between the calculated and measured values: 


Aa=a a. 


lob - lpom 

——tjb = 0.6 
——t/b = 0.8 
—retib = 1.0 


43 45 


o. [°] 


Fig. 22. Differences between measured and calculated exit angles a, vs. 
pitch t/b and stagger angle a, 


The diagram in Fig. 22 reveals that for large stagger angles 
in the N3 cascade, the differences between the measured and 
calculated angles are negligible. For smaller stagger angles, 
systematic increase of the difference between the calculated 
and real values is observed. Moreover, a remarkable effect of 
the cascade pitch on the range of these differences is observed. 
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VR GO 67 


0.95 1.0 1.05 1.1 1.15 12 


In extreme cases, for large t/b and small stagger angles the 
difference between the measured and calculated exit angles 
can approximately reach as much as 4.5°. 
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Operational causes of fatigue failures within 
passages of gas turbine engines 


Zbigniew Korczewski, Assoc. Prof. 
Gdansk University of Technology 


ABSTRACT 


In this paper a short presentation of fatigue wear process of structural elements of gas 

turbine engines has been given. The primary causes of fatigue crack formation within 

engine mechanical system and flow passages have been highlighted. Special attention 

was paid to low-cycle fatigue associated with unsteady heat-and-gas-flow processes 

developed in the passages. The selected damages have been demonstrated of gas flow 

paths of the engines operating in aviation, navy and power industry, along with origins 
of their formation and growth. 


Keywords: gas turbine engines, passages, fatigue wear, operational failures 


INTRODUCTION 


In the period of intensive investigations to improve 
operational reliability, durability and economic merits of 
gas turbine engines, a problem of effective methods of their 
technical state assessment becomes more and more important 
especially that in recent years a trend is observed of passing 
from the operation-based, planned maintenance to that based on 
current technical state, [1, 6, 7]. This is specially recommended 
when an unforeseen failure of engine is associated with a great 
hazard to traffic safety at sea, air and land [9, 14]. 

It is especially complicated to evaluate unserviceability 
states of turbine engine flow passages irrespective of an area 
of engine application. Technical state assessment traditionally 
performed on the basis of service measurements of gas-dynamic 
parameters does not bring satisfactory results in the case of 
fast developing surface defects of flow passages, which do 
not generate usually observable diagnostic symptoms. Service 
experience justifies that just most troublesome are flow passages 
of engine high-temperature part which is most exposed to 
thermal low-cycle fatigue whose consequences bring many 
hazards to engine safe operation. It is influenced by many 
factors, and the most important of them is associated with 
destructive and unavoidable effects of varying mechanical and 
thermal loads applied to elements of engines during their work 
in unsteady states [2, 3, 14]. 


LOW-CYCLE FATIGUE 


One of the parameters which impact turbine engine 
durability the most is fatigue strength of material of which 


its structure is built. During engine operation on a ship cyclic 
elastic and plastic deformations of structural elements take 
place, which result from multifold changes of mechanical loads 
(caused by gas forces generated by thermal dynamic cycle of 
engine as well as by rotary motion inertia forces of rotor units) 
and thermal loads (resulting from temperature gradients within 
structural elements). They cause variable internal stresses 
developing finally resulting in fatigue cracks. With a view of 
initial causes of fatigue failure development the following can 
be distinguished: 

- high-cycle fatigue (mechanical), which is characteristic for 
engine rotating units (engine rotor units) and occurs below 
yield point of structural material. Fatigue failures appear 
after a cycle number of load changes over 10* (usually over 
107) [3,4,5,9]; 

- low-cycle fatigue (thermal) which concerns the engine 
elements exposed to high-temperature effects due to exhaust 
gas flow (combustion chamber, stationary load-carrying 
elements of turbine part) which occurs over yield point of 
material. Fatigue cracks appear already below the cycle 
number of load changes equal to 10* [3, 4, 5, 10]. 


Relation between stress amplitude and cycle number of 
load changes, at which structural elements become failed due 
to material fatigue, is represented by the so - called Wóhler 
curve [5,10]. 

High-cycle fatigue of rotating elements of turbine engine 
rotor units is caused by vibrations, that means resonance 
phenomenon to appear: 

- in global sense when the entire engine suffers vibrations 
and then it constitutes a problem for designers, 
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- inlocal sense when one of the natural vibration frequencies 
of a given engine element, f, , appears close to the 
frequencies of periodically changeable vibration excitating 
forces, f. , QIn this case the increasing of deformations and 
stress amplitudes takes place and after a longer period of 
engine's operation in such conditions a crack may appear 
in the element and its failure may happen due to material 
fatigue. 


Maximum amplitude of resonance vibrations depends on 
stiffness and damping properties of supports of engine rotor units. 
The higher values of stiffness coefficients the smaller amplitude 
of resonance vibrations [1, 3]. 

Loss of stability of engine mechanical system under 
operation results from occurrence and development of the 
following detrimental phenomena: 

- Slowly worsening unbalance state of rotor units resulting from 
sediments, erosion, corrosion and bent shafts; 

- sudden increase of unbalance state of rotor units resulting 
from broken-off fragments of a rotating element, e.g. 
a blade; 

- exceedance of allowable load of bearings; 

- rotors seizing in engine body elements in labyrinth 
sealings, 

- appearance of self-exciting vibrations of rotor blades, 
mainly in compressors, (so-called flutter), during steady 
flow of working medium [11], 

- airandexhaust gas pressure pulsations within flow passages, 
as well as due to non-uniform distribution of temperature 
and flow velocity around flow passage circumference. 


Mechanical resonance phenomenon is of a very great 
influence on engine's serviceability and durability. However 
especially dangerous to its reliability is thermal fatigue of 
structural elements of combustion chamber and turbines, as 
a result of occurring unsteady processes such as start-ups, load 
changes (accelerations and decelerations), as well as engine 
turn-offs from operation. 

For instance, during engine's starting-up, within a few 
seconds structural elements of combustion chamber are exposed 
to dynamic thermal loads characterized by instantaneous 
exhaust gas temperature increments reaching about 80+90 
K/s [16]. In the case of exceeding the allowable temperature 
gradients, pumping phenomenon and even that of structural 
material creeping and deformation of turbine blades can occur, 
Fig. 1. 


Fig. 1. Deformations of LP turbine rotor blades, resulting from structural 
material creeping (acc. endoscopic examination) 
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In extreme cases, local burning-off edges and tips of HP 
turbine blades can happen, Fig. 2 and 3. Such failures of ship 
engine flow part usually result from sudden and unforeseeable 
operational events such as: 

* choking the exhaust gas outflow due to fire in outlet 
channel, 

* choking the engine air inlet (frozen air inflow shutters, 
sucking-in random things such as protective covers of 
shipboard devices etc). 


And, in the case of immediate stopping the engine during its 
operation in nominal (or near nominal) conditions, associated 
with eliminating idle running cooling phase, rotor seizure in 
turbine body may happen, Fig. 4 and 5. 


Fig. 2. Burnt-off tips of HP turbine rotor blades 
(acc. endoscopic examination) 


L 


Fig. 3. Burnt-off tips of HP turbine rotor blades 
(acc. endoscopic examination) 


Fig. 4. Increased tip clearance of HP turbine blade with visible traces 
of seizing in engine body (acc. endoscopic examination) 


direction v 
of observation 


“honeycomb“ 


Fig. 5. Wear of sealing of HP turbine (HPT) rotor blades due to their 
seizing in engine body (acc. endoscopic examination) 


It usually results in an intensive wear of body sealing 
of „honeycomb” type (Fig. 5), that in consequence leads to 
significant increasing the clearance of rotor blades of turbine 
and dropping its efficiency even by 10% [13]. It impairs efficient 
performance of entire engine because, as a result of increased 
radial clearance of HP turbine, character of distribution, to 
successive turbine stages, of disposed enthalpy drop becomes 
changed. The largest power drop occurs at the last turbine 
stages, i.e. engine propulsion turbine [12,13]. 

Another hazard to engine reliability, which results from 
rotor seizing in engine body, consists in possible cracking 
the compressor body, Fig. 6. The presented failure has been 
identified during endoscopic examination of a GTU6a single- 
shaft engine (of block structure) driving ship electric generating 
set. Compressor’s body was made of 40HNMA Cr-Ni low-alloy 
steel of a higher metallurgical purity. The situation occurred as 
a result of immediate switching-off the engine from operation 
under nominal load and passing over the cooling phase under 
idle running. Such decision of operators resulted from a fire 
observed in exhaust gas channel and entering the engine into 
unsteady operation range. After that the engine was dismounted 
and sent to its producer. Repair of the body consisted in vacuum 
electron-beam welding the crack. 

On the basis of the character of the appeared failure of the 
compressor body it can be concluded that the initial cause of the 
crack was a material defect, e.g. presence of sulfides at place 
of fracture. In effect of action of dynamic loads the fatigue 
fracture propagating just from the edge of the material defect, 
was initiated. As a result of rotor seizing in the compressor body 
in the temperature of about 700-800 K, further propagation of 
the gap took place in the form of cleavage fracture. As the body 
was cooled more fast than the rotor the gap became separated 
mainly due to tensile forces (axial ones) with an influence of 
shear forces. 


Designers of gas turbine engines take into account possible 
stopping the engine working under high load with passing 
over the cooling phase under idle-running. It is covered by 
manufacturer’s testing program. Therefore the seizing of rotor 
units in engine body should not happen at all. However in the 
described case the engine immediate stopping was preceded by 
its unsteady work when a flash-back towards engine intake took 
place causing an inadmissible, very intensive rise of compressor 
material temperature much over its permissible value of 520 K, 
determined by the manufacturer [15]. 


Fig. 6. Internal surface of GTU6a turbine engine compressor body, showing 
a fatigue fracture (acc. endoscopic examination) 


The presented examples, despite their immediate destructive 
consequences to the engine’s reliability, had also secondary 
consequences to its durability. During every load change, 
deformations and thermal stress changes of engine structural 
elements took place. In elevated temperatures they exceed 
material yield point and cause plastic compression-tension 
deformations at every heating-cooling cycle of the material. 
It turns out that even after a small number of such cycles the 
material surface cracking of deformed elements can happen, 
which, while propagating into their structure, can lead to 
dangerous breaking-off a fragment of the structure. This is an 
especially hazardous defect in the case of combustion chamber 
flame tubes, Fig. 7. There are known initial causes of an aircraft 
crash of tragic consequences where due to low-cycle fatigue of 
flame tube material of one of the engines fatigue cracks in its 
structure occurred. During start-up of the engine a relatively 
large metallic fragment was broken out of it. The fragment 
protruded engine casing and wing sheeting and caused fuel to 
leak and in consequence a violent fire to break in the rear part 
of the fuselage, that took the toll of 60 persons [acc. Discovery 
Science]. 

If low-cycle fatigue characteristics of a structural material 
are known it can be possible to easily determine a number of 
cycles after which surface cracking may happen, depending on 
a single amplitude value of the strain e which occurs during one 
heating-cooling cycle of the material operation, Fig. 8. 

From the numerical data given by the fatigue characteristics 
the conclusion can be made that if the engine, for example an 
UGT3000 Zorya engine, has to be started 2+3 thou times during 
its life-time determined to be 12 years equivalent to 4 thou 
working hours, then one-sided deformations of the material 
(short-range local strain) should not exceed 0.2+0.3 %, at each 
start-up, [15]. 

Stress concentration occurring in places of significant 
curvature changes of structural element outer surfaces (i.e. 
at notches) is conductive to local strains and fatigue fracture 
generation. Hence they appear mainly in necks of load-carrying 
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Fig. 7. Low-cycle fatigue fracture of flame tube of turbine engine 
combustion chamber: a) picture recorded just after the engine S 
dismounting; b) picture obtained by endoscopic examination 
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Fig. 8. Low-cycle fatigue characteristics: € — strain (single amplitude), 
N — number of cycles (heating - cooling), a, b — fatigue limit 
(shown for two different materials) 


disks of rotors, in places of the greatest diameter changes, close 
to openings and gaps of combustion chamber flame tubes as 
well as in neighborhood of significant changes of profiled 
cross-sections of internal shell bodies, Fig. 9. 
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Fig. 9. Low-cycle fatigue fracture on internal body surface, under HP. 

turbine guide vanes of a Zorya UGT3000 engine: a) picture recorded 

just after the engine s dismounting; b) picture obtained by endoscopic 
examination. 


CONCLUSIONS 


Crucial diagnostic problem of gas turbine engines is to be 
able to early detect and correctly identify initial failures of 
structural elements of engine flow part during early phase 
of their development. 


While analyzing the research results presented in this 
paper it seems impossible to overestimate importance of 
new possibilities introduced, to the aspect in question, 
by more and more dynamically developing endoscopic 
methods. 


On the basis of available statistical data as well as results 
of the author's research it can be concluded that in present 
about 60% of all recorded failures of flow part of turbine 
engines is detected by means of endoscopic methods, which 
are also identified by other diagnostic methods. 


However in the case of fatigue failures endoscopy 
effectiveness is even higher, provided the diagnosing 
personnel is well experienced and possesses deep knowledge 
on probable places of occurrence and characteristic features 
of surface defects which have been so far happened in 
engines of given type series. 


BIBLIOGRAPHY 


. Bolinski B., Stelmaszezyk Z.: Aircraft drives. Operation 

of turbine engines (in Polish). WKiL (Traffic and 
Telecommunication Publishers), Warsaw 1981 

. Cohen H., Rogers G.F.C., Saravanamuttor H.I.H.: Gas turbine 
theory. Longman Scientific & Technical, New York 1987 

. Dzierzanowski P. et al.: Aircraft drives. Turbine propeller 

and helicopter engines (in Polish). WKiL, (Traffic and 
Telecommunication Publishers), Warsaw 1985 

. Kocanda S., Szala J.: Background of fatigue calculations (in 
Polish). PWN(State Scientific Publishing House), Warsaw 1985 
. Kocanda S.: Fatigue fracture of metals (in Polish). WNT 
(Scientific Engineering Publishers), Warsaw 1985 

. Korczewski Z.: A method for diagnosing flow part of ship gas 
turbine engine in service (in Polish). Doctor thesis. AMW(Polish 
Naval Academy), Gdynia 1992 

. Korezewski Z.: Diagnostic-purpose identification of gas- 
dynamic processes in gas turbine engine compressor unit (in 
Polish). AMW (Polish Naval Academy), Gdynia 1999 

. Korezewski Z.: Endoscopy of ship engines (in Polish). AMW 
(Polish Naval Academy), Gdynia 2008 

. Lewitowicz J., Wieczorek O. Zyluk A.: Damages of structural 
components and units of a gas turbine engine. Journal of Polish 
CIMAC, Diagnosis, Reliability and Safety, Vol. 4, No. 2/ 2009 


10.Neimitz A.: Fracture mechanics (in Polish). PWN(State 


Scientific Publishing House), Warsaw 1998 


11.Rzadkowski R.: Flutter of turbine rotor blades in inviscid flow. 
AMW (Polish Naval Academy), Gdynia 2004 

12.Szezecinski S.: Aircraft two-rotor and two-flow turbine engines 
(in Polish). WKiŁ (Traffic and Telecommunication Publishers), 
Warsaw 1971 

13.Szezecinski S.: Study on tip clearance of rotor units of aircraft 
turbine engines, taken as a design and operation parameter 
(in Polish). Appendix to Bulletin of WAT(Military Technical 
Academy), No.4 (248), Warsaw 1973 

14.Soares C.: Gas turbines. A handbook of air, land and sea 
application. Butterworth-Heinemann, Elsevier Inc. 2008 

15.Technical and operational documentation of ship gas turbine 
engines of the types: GTU6a, DE59, UGT Zorya, LM2500 
General Electric. 

16.Reports on diagnostic tests of piston and turbine combustion 
engines used on Polish Navy ships (in Polish). Research 
projects, AMW (Polish Naval Academy), Gdynia 1992 + 2008. 


CONTACT WITH THE AUTHOR 
Zbigniew Korczewski, Assoc. Prof. 
Faculty of Ocean Engineering 
and Ship Technology 
Gdansk University of Technology 
Narutowicza 11/12 
80-233 Gdansk, POLAND 
e-mail : sek4oce@pg.gda.pl 


POLISH MARITIME RESEARCH, No 1/2010 61 


POLISH MARITIME RESEARCH 1(63) 2010 Vol 17; pp. 62-68 
10.2478/v10012-010-0007-2 


Diagnostics, maintenance and regeneration 
of torsional vibration dampers for crankshafts 
of ship diesel engines 


Wojciech Homik, Ph. D. 
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ABSTRACT 


Periodically changeable gas and inertia forces which occur during operation of engine 
generate transverse, axial and torsional vibrations of crankshafts of multi-cylinder 
combustion engines. Torsional vibrations are those which endanger crankshafts of 
multicylinder combustion engines the most. In order to minimize their impact a torsional 
vibration damper is installed at crankshaft s free end. Its technical state directly influences 
lifetime and reliability of engine. In this paper methods of diagnosing, maintenance and 
regeneration of torsional vibration dampers used in shipbuilding, are discussed. Also, are 


presented results of multi-year statistical investigations carried out in cooperation with a firm maintaining 
and regenerating ship engine torsional vibration dampers, which illustrate types of failures occurring in 
viscous and spring torsional vibration dampers. 


Keywords: crankshaft, torsional vibrations, damping of torsional vibrations, torsional vibration dampers, 
durability and reliability of multi-cylinder engine, diagnostics, maintenance, regeneration 


INTRODUCTION 


Multi-cylinder combustion engine in operation generates 
vibrations which directly and detrimentally impact lifetime of 
engine’s parts and the whole propulsion system. 

Crankshaft is one of many engine’s parts which is exposed 
to effects of forced vibrations. Periodically changeable gas 
and inertia forces occurring in piston combustion engines 
generate: 

- transverse vibrations (Fig. la) 

- axial vibrations (Fig. 1b) 

- torsional vibrations (rys. 1c) 

which produce, by affecting flexible engine crankshaft, its 
changeable deformations [1, 2]. 


Transverse vibrations produce bending deformations of 
crankshaft between its support points, i.e. main bearings. 
However they do not endanger the cranklshaft seriously 
because it has a low transverse flexibility resulting from the 
way of its bearing in engine’s body. Practically every crank of 
the shaft is supported by main bearings. Due to crankshaft’s 
design solution itself also axial vibrations affect very seldom its 
operation especially in the case of engines used for propulsion 
of land vehicles. However such vibrations are able to generate 
failures of ship propulsion engine more often. 

Resonant torsional vibrations are most dangerous for 
engine’s crankshaft [1,2,3], as - in contrast to transverse and 
axial vibrations - they do not propagate to other parts of the 
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Fig. 1. Kinds of crankshaft vibrations: a) transverse, b) axial, c) torsional 


engine, e.g. bearing casings, and in many cases they also do not 
generate noise which is a factor informing engine’s operator on 
an incorrect work of the engine. It results from that torsional 
deformations of cranks of the shaft are generally not restrained 


by any factor except its rigidity hence torsional vibration 
amplitude can exceed allowable limits leading to twist damage 
of the crankshaft. 

Users of combustion engines expect the engines used 
by them to be reliable and durable. Engine designers make 
all efforts to place frequency of k-th harmonic of excitation 
force as far away as possible from natural torsional frequency 
of engine’s crankshaft. In ship large multi-cylinder engines 
the problem can not be fully solved. For design and service 
reasons changes of engine rotational speed (service speed), run 
of excitation forces, or natural frequency of shaft are difficult 
to be made [4]. In such cases combustion engine designers 
decide to introduce an intermediate solution and apply vibration 
eliminators, i.e. torsional vibration dampers which are installed 
at free end of engine crankshaft (Fig. 2). 


uL. ;- EIN 


Fig. 2. a) 8S70MC-C diesel engine for ship propulsion [18], 
b) Crankshaft of RTA96 Wartsila-Sulzer diesel engine, 
the largest-in-the- world [16], 1- torsional vibration damper 


Their task is to lower amplitudes of torsional vibrations of 
engine crankshaft. A correctly designed damper of torsional 
vibrations makes it possible to lower even ten times resonant 
amplitude of torsional vibrations [5,6,7]. However it should 
be remembered that every damper consumes a part of engine's 
effective output. 

For years to minimize hazard resulting from torsional 
vibrations in ship diesel engines the following types of dampers 
have been applied: 

- frictional dampers, 
- rubber dampers, 

- viscous dampers, 

- spring dampers. 


The above specified dampers are typical dynamic ones in 
which inertia force is used to damp torsional vibrations. 

Despite their common name the dynamic dampers differ 
to each other not only by design solutions but also operational 
characteristics (Fig. 3) [8]. 

Irrespective of an applied design solution to dynamic 
damper, its technical state is decisive of its effectiveness, which 
greatly depends on service conditions of engine at which the 
damper is installed. As majority of torsional vibration dampers 
operate in periodically changeable conditions (e.g. changeable 
rotational speed) or changeable atmospheric conditions the 
dampers should be periodically diagnosed like other devices 
to ensure correct operation of engine. 

Practically it is not possible to elaborate universal 
diagnosing methods for all kinds of dampers of ship engine 
torsional vibrations. This mainly results from differences in 
their design. 

Every producer of dampers should have at his disposal 
effective methods for diagnosing vibration dampers to make 
it possible to check if a given used damper still maintains its 
effectiveness. One of the damper’s effectiveness criteria can be 
magnitude of torsional vibration amplitude of the shaft at which 
the damper is installed: |E max| € € dop. However the criterion can 
be insuffcient for assessing damper's effectiveness. Taking the 
fact into consideration, the damper producers should require 
- from the side of operators of the engines at which torsional 
vibration dampers are installed - to perform periodical diagnostic 
tests and overhaul of the dampers with help of personnel of 
specialist firms or certified servicemen. Irrespective of a design 
solution, all torsional vibration dampers fixed on ship engine 
crankshafts should be subjected to periodical survey every 
18000 + 20000 hours of operation. 


DIAGNOSTICS OF FRICTION TORSIONAL 
VIBRATION DAMPERS 


Friction torsional vibration dampers of Lanchester type 
have been the first applied to damping torsional vibrations of 
ship engine crankshafts (Fig. 4). 

In such dampers the inertia ring's motion relative to the 
boss rigidly connected with the engine crankshaft is possible 
to be realized only when the inertia force moment M, of 
the active inertia I, is greater than the friction moment M.. 
During the motion resulting from dry friction between friction 
linings, torsional vibration energy is transformed into heat and 
discharged to the environment [9]. 

Motion of the inertia ring relative to the fixed boss, due to 
high frequency of vibrations and the amplitude of crankshaft 
torsional vibrations Q.... > Pi results in an intensive wear of 
friction lining which leads to changing value of friction force 
and the so - called detuning of the damper. 
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64 


a) M 


b) 


Fig. 3. Operational characteristics of dampers: a) Lanchester 5 frictional damper, b) a damper of non-linear characteristics, 
c) a spring dynamic damper of linear characteristics, d) a spring damper of prestrained springs 
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Fig. 4. Friction torsional vibration dampers of Lanchester type: 
1) shaft, 2) boss, 3) inertia ring, 4) disc, 5) friction lining, 


6) tension spring, 7) driving dog 
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The fact directly impacts effectiveness of the damper and 
forces the engine operator to perform its frequent monitoring 
and maintaining. 

Control of friction torsional vibration damper is limited 
to checking quality state of friction surfaces and regulation 
of holding down force of linings with the use of the tension 
springs. In justified cases, i.e. when wear of friction lining 
exceeds a limiting value, producers ofthe dampers recommend 
definitely their replacement. 

Lack of possibility to maintain value of the friction 
moment M, constant during operation of the damper results 
from difference in rate of wear of friction linings as well as 
from getting water, oil or mechanical contaminations between 
friction surfaces. Despite their doubtless merits the friction 
torsional vibration dampers have fallen short of expectations 
put in them and have been soon replaced by other dampers. 


DIAGNOSTICS OF RUBBER TORSIONAL 
VIBRATION DAMPERS 


In shipbuilding industry rubber torsional vibration dampers 
were first time used in 1915. Big hopes have been put in them 
but service experience has verified it within short time. 

It has turned out that structure of rubber undergoes ageing, 
like many other organic materials, which results in increased 
viscosity, hardness or brittleness [10, 11]. Rubber ageing process 
results from action of oxygen as well as ozone contained in air 
and it very fast develops in high temperature. As a result of 


ageing rather large decrease of rubber strength properties as 
well as significant drop of its resistance to periodically changing 
loads takes place. 

Lack of stable physical properties of rubber, connected with 
the above described phenomena, has not provided the rubber 
torsional vibration dampers with a greater role in shipbuilding 
industry. At the end of the 1950s application of rubber dampers 
was practically terminated. In spite of that they have been still 
in use on some older floating units. 

Their short lifetime does not up to now prevent designers of 
car engines from their successful application. It mainly results 
from that car engine’s lifetime is much lower than that of ship 
engine and is equal at the most to about 2000000 km (trucks, 
buses). As results from tractional tests of car combustion 
engines fitted with rubber torsional vibration dampers their 
operational effectiveness drops after about 180000 km of 
mileage. 

The rubber torsional vibration damper consists of the boss 
and plunger connected to each other with the use of the rubber 
ring (Fig. 5). The rubber ring has an appropriate hardness, 
flexibility and internal damping. 


Fig. 5. A rubber torsional vibration damper undergoing overhaul in 
DAMPOL Co.: 1) inertia ring, 2) boss, 3) rubber damping element 


Irrespective of damper’s design solution and application 
(either to ship or car engines) their producers recommend to 
perform periodical survey mainly from the point of view of 
technical state of rubber. Producers of rubber dampers consider 
them to be disposable, i.e. those which should be replaced 
by new ones after a given period of service (ship engines) or 
mileage (car engines). As a rule their producers do not provide 
for any repair or regeneration of such dampers. 


DIAGNOSTICS AND REGENERATION 
OF VISCOUS TORSIONAL VIBRATION 
DAMPERS 


In 1929 in U.S. shipbuilding industry viscous torsional 
vibration dampers were applied to damping torsional vibrations 
of governors of 3000 HP engines installed on submarines 
[12]. However no satisfactory results were obtained [13] 
though an innovative, as in those days, solution consisting 
in parallel connection of two dampers filled with different 
viscosity liquids, was used. The main reason of the failing 
was application, as damping medium, of silicone oil produced 
from organic compounds, whose viscosity was dramatically 
decreasing along with the increasing of operational temperature 
of the damper. The high operational temperature had also 
decisive impact on oil ageing rate, hence also on the damper’s 
reliability and lifetime. The doubtless drawback of the viscous 
torsional vibration dampers resulted in that for many years they 
have been not applied at all. 

They experienced their renaissance at the end of the 1950s 
when Dow Corning Corporation introduced engine oils on the 
market, whose physical properties ensured to obtain appropriate 
technical parameters of dampers in question. The improved 
effectiveness of the ,,new" damper made that it squeezed the 
friction and rubber dampers, out of the market, which have been 
applied so far, especially in shipbuilding industry. 

The viscous torsional vibration dampers are composed of 
the tight casing filled with silicone oil, in which the inertia ring 
floats (Fig. 6) [5]. 


a view on a marine viscous torsional vibration damper [17]: 1) casing, 
2) axial bearing, 3) inertia ring, 4) plug, 5) journal bearing, 6) oil 
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Diagnostics of such dampers depends on their design 
solutions. In present two design solutions are applied in 
industry: non-dismountable dampers (in car industry) and 
dismountable dampers (in shipbuilding industry). 

Diagnostics of viscous torsional vibration dampers consists 
first of all in checking the vibration amplitude criterion and then 
measuring value of damping liquid viscosity [13]. 

Testing samples of the liquid (oil) should be taken from 
the damper being in neutral environment in order not to cause 
connection of oxygen and hydrogen atoms and tearing this way 
particle bonds and lowering real viscosity of the oil [14]. In 
the case of stating an inappropriate viscosity of the tested oil 
its replacement is recommended after very accurate cleansing 
all parts of the damper. 

Statistical tests showed that the decreasing of oil viscosity 
during operation of dampers is a natural factor (Fig. 7). In many 
cases during diagnosing it turned out that it was not possible 
to take oil samples as its actual viscosity was several times or 
even a few dozen times greater than its initial viscosity (Fig. 8). 
The so significant increase of oil viscosity leads to the inertia 
ring immobilizing, hence to the situation in which the damper 
becomes a vibration exciter. 
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Fig. 7. Statistical results of testing the dampers fitted with bronze bearings, 
v [?6] = V oa imitar 100% 
Kind of material used for inertia ring and its journal 
bearings and axial ones greatly influences changes of silicone 
oil viscosity and thereby damper’s lifetime [14]. 


Fig. 8. Used silicone oil (solidified material) stuck to the casing of the 
viscous torsional vibration damper 
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Regeneration of viscous torsional vibration dampers 
consists in a.o. replacing their journal and axial bearings made 
of bronze with the bearings made of PTFE teflon. 

Effectiveness of such operations is confirmed by results of 
statistical tests shown in the diagrams of Fig. 9, from which it 
results that in the case of replacing bronze bearings with teflon 
ones the above described phenomenon of the large increase of 
oil viscosity due to diffusion occurring at the place of contact 
between bearing and casing, was almost entirely eliminated, 
thereby mechanical failures of inner surfaces of the damper 
were prevented. 
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Fig. 9. Statistical results of testing the dampers fitted with teflon bearings, 
v[?6] =v /V 
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Fig. 10. Failures of inner surfaces of a viscous torsional vibration damper 


In the case if mechanical failures on inner surfaces are 
disclosed during diagnosing the damper (Fig. 10) their 
regeneration is strongly recommended (turning the damaged 
surfaces and their repair by welding in justified cases). During 


regeneration of inner surfaces of a damper, changes of their 
geometrical parameters occur thereby also changes take place 
of dynamic properties of the damper. In such cases the damper 
should be filled with a liquid of a greater viscosity than that initial 
in order to restore appropriate effectiveness of the damper. 

Tightness of the damper consitutes a factor of great 
influence on its correct operation. During overhaul of the 
damper its sealing is replaced and after its closing and filling 
its tightness test is performed. 

Also, non-dismountable viscous torsional vibration 
dampers, like the above described rubber dampers, are not 
subjected to repair and regeneration and considered by their 
producers to be disposable parts. 


DIAGNOSTICS AND REGENERATION 
OF SPRING TORSIONAL VIBRATION 
DAMPERS 


Research on a spring torsional vibration damper was initiated 
in 1958 when it turned out that viscous vibration dampers did 
not effectively lower crankshaft torsional vibrations in power 
system engines. The first spring damper was manufactured 
by Geislinger firm and installed on the Danish ferry ,,Holger 
Danske” in 1962. The spring torsional vibration damper, like 
rubber one, is composed of a boss and inertia ring mutually 
connected with the use of flexible elements in the form of spiral 
springs, packets of flexible pads and bush springs. 

Today such dampers can be commonly found in ship 
propulsion engines as well as large power plants. 

Two types of spring dampers are presently produced for 
shipbuilding industry: 
- Bush spring damper acc. Pielstick design (Fig. 11) 
- Damper with packets of straight springs (radial springs) 

acc. Geislinger design (Fig. 12). 


As compared with viscous trosional vibration dampers the 

spring dampers are characteristic of: 

- smaller dimensions, 

- small inertia moment of inertia ring relative to mass inertia 
moment of boss, 

- greater resistance to mechanical failures, 

- greater lifetime, 

- higher allowable temperature of operation. 


Fig. 11. Ship engines spring torsional vibration damper with double packet 
of torsional springs, of Pielstick design: a) general view, b) socket with two 
packets of springs, c) socket with single packet of springs, 

d) packet of springs 


Spring dampers should be serviced once a year at least 
and their overhaul should be performed after about 20000 h 
of operation. Apart from checking the basic effectiveness 
criterion (amplitude criterion) of such damper its overhaul 
consists mainly in: 

- control of its tightness (after about 2000 h of operation), 

- control of patency of lubricating channels through which 
oil is delivered between cooperating surfaces of boss and 
inertia ring, 

- examination of technical state of flexible elements (springs, 
packets of springs) from the point of view of possible 
mechanical failures, 

- control of technical state of bronze intermediate sliding 
pads. 


In spring dampers, depending on design solutions, 
geometrical parameters of the sockets in which springs or 
packets of springs are placed, should be also controlled, 
especially in the case of the dampers with prestrained springs. 
In the case of the dampers with packets of radial spring plates 
(Fig. 12) the side clearance ,,Z" between side surfaces of springs 
and sockets should be checked (Fig. 13). 


i 


Fig. 12. Spring torsional vibration damper with packet 
of radial springs, of Geislinger design 


Fig. 13. Side clearance ,,z” 


Values ofthe clearances depend on damper's size and reach 
a few dozen micrometers. 

Statistical tests confirm that spring torsional vibration 
dampers are characterized by high reliability and their 
diagnostics consits mainly in adjustment of magnitude of 
clearances. (Fig. 14). 
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Fig. 14. Statistical results of testing the spring torsional vibration dampers 
CONCLUSIONS 


e Repair or regeneration of all the above discussed dampers 
is associated with their disassembling off the crankshaft at 
which it is fixed, and their dismounting into parts. Before 
starting the damper’s disassembling its position relative to 
the crankshaft should be marked and before its dismounting 
mutual position of all its parts should be also appropriately 
marked. 


e Precise technical parameters dealing with diagnosing 
and servicing the torsional vibration dampers are kept 
confidential. They are only at disposal of producers of 
the dampers as well as firms which are engaged in their 
maintenance and regeneration under supervision of 
the worldwide ship classification institutions such as: 
Germanischer Lloyd, det Norske Veritas, Lloyd’s Register, 
Bureau Veritas, Polish Register of Shipping, Russian 
Register of Shipping. 


e The presented results of testing the torsional vibration 
dampers carried out in cooperation with their producer 
[15] prove that majority of users of combustion engines 
fitted with torsional vibration dampers do not appropriately 
monitor their technical state. It leads to the situation in 
which a vibration damper is converted into a generator of 
vibrations or undergoes damage (Fig. 15). 


Fig. 15. A damaged spring torsional vibration damper of Pielstick [15] 
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Thorough capital repair of torsional vibration dampers are 
mainly performed during capital repair of entire ship propulsion 
systems. 
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Correction of the combustion engine run 
irregularity in hybrid systems 
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ABSTRACT 


Possibilities ofcombustion engine — electric machine assembly run irregularity improvement 
have been presented in the following paper. Such system is the basic one in numerous hybrid 
car available on the market but the concept of using an electrical machine for reduction 
of rotational speed changeability amplitude in reciprocating machines could be utilized in 
any drive unit, especially in ship propulsion. This idea consists in generation of braking 
torque by the electric machine whenever the coupling shaft increases its rotational speed 
and vice versa, driving the system when the momentary speed is lower than the average 
speed within a complete shaft rotation. 


Keywords: hybrid drive, run irregularity 


INTRODUCTION 


From the very beginning of its history an internal 
combustion engine which is one of the most important 
inventions initializing technical civilization created problems 
resulting from the changeability of torque transferred to 
the power receiver. Contrary to electric motors or turbines 
cyclicity of thermodynamic changes which cause the cyclicity 
of generated torque is the essence of reciprocating machine 
operation. However, there is no thermal machine as efficient 
as the piston engine in energy transformation from chemical 
to mechanical one. A question arises whether the basic 
disadvantage of reciprocating engine consisting in changeability 
of instantaneous speed within a single cycle could be eliminated. 
A flywheel, invention known from the very beginning of engine 
history serves as the device moderating the speed irregularity. 
Certainly, an increase in cylinder number contributes to the 
speed equalization as well. Unfortunately, both solutions used in 
almost every combustion engine bring about specific drawbacks. 
The flywheel increases weight of engine and vehicle as well. 
As a matter of fact a similar effect can be achieved by the 
increase in flywheel overall dimensions but hardly ever engine 
dimensions are defined by the flywheel. As the result dimensions 
and moment of inertia are the compromise of the possibly lowest 
speed irregularity and engine weight and dimensions. Without 
going into details one can assume that for a defined engine swept 
volume an increase in cylinder number leads to the deterioration 
of mechanical efficiency and to an increased specific fuel 
consumption as well. Other group of problems connecting with 
the subject of engine run irregularity are the torsional vibrations 


within the range of unit natural frequencies. Following the 
development of combustion engines one can notice that the 
problems of resonance vibrations were the prime mover of the 
research on limitation of instantaneous speed change amplitude. 
Insufficient recognition of resonance vibration problem usually 
ends in a user dangerous failure consisting in the shaft breakup. 
The only solution that allow to avoid hazards is the use of 
torsional vibration dampers. Unfortunately, damping means 
energy dissipation and unavoidable increase in fuel specific 
consumption. 

So what can be done to avoid the above mentioned 
disadvantages of solutions aimed at improvement of engine 
run irregularity? 

The answer could be an application of hybrid drive 
to automobiles. There are forecasts that in a near future 
a considerable number of cars will be equipped with an 
arrangement facilitating the recovery of energy during braking 
and then using it to accelerate the car. However, the concept of 
application of electrical machine to the reduction in rotational 
velocity of reciprocating machines could be utilized in any drive 
unit including ship propulsion in particular. Such arrangement 
requires use of a set consisting of: 


* generator that can operate as electrical motor, 
* converter modifying the generated current to battery, 
* suitable batteries. 


In order to use the arrangement reasonably its power 


should correspond to the power of combustion engine, at 
least it should not be less than 1096 of engine nominal power. 
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Tab. 1. Parameters of electrical energy accumulators where energy is freed as a result of reversible chemical reactions 


Power 
concentration 


Energy 
concentration 


Battery type 


kWh/kg 
0.03-0.04 


Lead 0.2-0.3 


Life span 


Cost of energy 
concentration 
unit 


Time of energy 
and power 
preservation k 


Number of charge 
cycles 


zi per kWh 


300-400 400-600 


Nickel-Cadmium 0.04-0.05 0.080-0.180 


>2000') 2 500 


Nickel-Metal 
Hydride 
(Ni-MH) 


0.06-0.15 0.20-0.30 


>1000 1200 


Sodium-Nickel- 
Chloride 
(NaNiCl) 


0.09-0.10 


Lithium-Ion 


(Li-Ion) 0.09-0.14 


0.30-0.60 


500-750 1000-2000 


Lithium-Polymer 


(LiPo) 0.11-0.13 


500 Data unavailable 1000 


Zinc-Air 0.10-0.22 


Data unavailable Data unavailable 250 


0.005 
0.20 


Capacitor 


Target values 0.50 


500 000 10 
1000 10 


Very high 
100 


8.5 
or about 
2?) 


Gasoline tank limited by engine 


power 


2 
or about 
10°) 


Practically no 


No limits e 
limits 


1) on condition that user prevents so called memory effect caused by charging battery not fully discharged; 
?) converted to the mechanical energy unit, which results from assumed efficiency of IC engine of about 2596 


Amount of energy gathered cyclically in battery is not great. 
A simple calculation proves that amount of energy required for 
acceleration of a passenger car is far lower than that contained 
in start battery, 1.e. about 1 kWh. As it turns out, so called 
energy concentration of known batteries is a crucial limitation. 
A classic lead battery of 100 Ah capacity and 12 V voltage can 
give back energy with satisfactory efficiency with power not 
exceeding 1 kW. So the use of batteries of higher capacity, e.g. 
1000 Ah becomes a necessity in the hybrid arrangements. On 
that score a good solution is the use of high energy capacitors. 
High power concentration and low energy concentration are 
their properties. Table 1 presents basic properties of electrical 
energy accumulators, including high energy capacitors [1, 2]. 

The idea of hereby study is an evaluation of possibilities 
to make use of car hybrid drive in order to improve the engine 
run irregularity. In the case of ship propulsion the control of 
torque variability of the auxiliary electrical machine coupled 
to the piston engine of main propulsion in order to remove 
the rotational speed variable component is far easier thanks to 
the lower frequency of this component in comparison to the 
automobile drives. 


AN INSTANTANEOUS POWER 
OF IC ENGINE 


A piston combustion engine generates momentary power 
far higher than the average one, i.e. power given by producer as 
nominal power. In principle, by the increase in cylinder number 
of an in-line engine one can achieve the mean power close to the 
instantaneous maximum power. Such direction of modifications 
leads to the improvement in the 6 engine run irregularity. Alas, 
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contrary situations also happen. Fig. 1 presents a course of 
torque momentary value generated by a 4-cylinder engine, 
while Fig. 2 presents the same for a 5-cylinder engine. 


Zbior — 35R410 

Le=4  om-367 Dt=0.081 r=0.048 etam=0.900 

eps=18.50 fic=1.73 ro=2.00 del=1/127 Msr=237 pa=0.093 
mp=1.01 mo=2.00 kz=0.200 lam=0.297 mi=0.000 id-0.107 


Moc Ni = il Max. norm. — 0.816 DI 
Moc/V - 53.04 [KA] Max. pred. — 18.261 [m/s 
Max. nor. — 4.207 [kN] czas obrotu — 17.1 [ms] 
Temp. sil. — 100.0 [C] 
— omega [rad/s] 
— moment [Nm] 
Torque [Nm] omega [rad/s] 


180.0 
0.0 8.6 17.1 25.7 


360.0 540.0 Crank angle. deg 
Time [ms] 


Fig. 1. Course of torque generated by an in-line 4-cylinder engine — green 
line and the course of speed produced by this torque 


The run irregularity level 6 marked as “del” and the mean 
indicated torque marked as “Msr” can be found above charts in 
3" row. The relation of torque maximum value to its mean value 
is 2.4 for a 4-cylinder engine and 3.5 — for the 5-cylinder one. 


Zbior — 35R510 
Le=5 om=367 Dt=0.081 r=0.048 etam=0.900 
eps=18.50 fic=1.73 ro=2.00 del=1/79 Msr=296 pa=0.093 


mp=1.01 moz2.00 kz=0.200 lom=0.297 mi=0.000 pd=0.107 
Moc Ni — 108.4 [kw] Max. norm. — 0.816 A 
Moc/V — 53.94 [KM/] Max. pred. — 18.261 [m/s 
Max. nor. — 4,207 (5 ] cezas obrotu — 17.1 [ms] 
Temp. sil. — 100.0 [C] 
— omega [rad/s] 
— moment [Nm] 
Torque [Nm] omega [rad/s] 
900 22 
700 18 
500 14 
300 10 
100 6 
eo S PENES D» GENNA GEO CD GERGEND QU Li NESS 1! ZEN 
"h o jis VINI 7 cm LI LV Mm mu Sn, 


WEN." SIAE E e RIO ION SC Go. 000 O NE Li 
PAR. PAL. MD a AD FINENDO LCUND Mic. CANEED e MATS 


0.0 180.0 360.0 540.0 Crank angle, deg 
0.0 8.6 17.1 25.7 Time [ms] 


Fig. 2. Course of torque generated by an in-line 5-cylinder engine — green 
line and the course of speed produced by this torque 


Presented examples show complexity of conditions affecting 
torque transmitted to the power receiver. Unfortunately, 
changeability of torque brings about negative results in form of 
vibrations which significantly increase the stress in power drive, 
i.e. in couplings, joints, gears and so on. These phenomena have 
been known for years. A counteraction consists in definition of 
vibration amplitude at resonance caused by a certain harmonics 
of torque and eventual check if admissible stress has not been 
exceeded in parts of drive train. Numerical methods used 
nowadays allow to predict parameters of vibrations generated 
simultaneously by all harmonics at any moment of engine 
run. Examples of such considerations will be presented in the 
next chapter. 


PARAMETERS OF ENGINE TO POWER 
RECEIVER COUPLING 


Moment of force and instantaneous speed are the principal 
parameters of engine to power receiver coupling. On a real 
object these parameters vary independently of any rules. 
Hence, classic method of determination of natural frequencies 
produced by a specific harmonics does not allow to define 
the probability of exceeding permissible parameters with 
satisfying accuracy. The process of vibrations develops in time 
and only determination of real parameters of the phenomenon 
in sufficiently long time secures unfailing transmission of the 
torque to a power receiver. Therefore exists a need for definition 
of torque and instantaneous rotational speed in order to predict 
possible problems in the drive train. Definition of parameters 
of drive transmission using the numerical methods seems to be 
a simple task. Fig. 1 presents a schematic of drive train from 
IC engine to power receiver. 

The mass moments of drive train elements inertia are 
marked as ©, the reduced length of elastic elements connecting 
inertial elements — L while the dumping between inertia 
elements — C. The key elements of drive train are: 

e crankshaft of the ®,, L,, and C,,, i.e. reduced moment, 
reduced length and damping, respectively, 

e flywheel of inertia moment ©, — stiffly coupled to the 
crankshaft, reduced length L,, and dumping C,,, and 
transmitting the torque to an additional flywheel of the ©, 
inertia moment, and so on. 


e C333 C34 
C, 

e 

e e e 
e uL. & 
m| 123 
= 33 
|a L34 = 


Fig. 3. Schematic of torque transmission 
from combustion engine to power receiver 


The interrelations facilitating determination of coupling 
torque and momentary speed of drive train elements with an 
assumed time step de/t are presented in procedure | [3]. 

In procedure 1: 

e variable b/j,i/ is the coordinate of individual mass revolution 
angle, 

* variable e/j/ is the yield resulting from shaft length Lili 
— see Fig. 3, 

e variable fer[j] is the mass moment of inertia ©, 

* flumj — dumping Cj,j+l, 

* whereas the sign of deltmo moment decides if moment is 
transmitted to the electrical motor which happens during 
stroke of expansion, or it is transmitted to the combustion 
engine, when auxiliary strokes are performed. 


Other terms of the procedure are not important to understand 
the method of correction of the power transmission shaft 
instantaneous speed. 


CONCEPT OF VIBRATION LIMITATION 
USING THE HYBRID ARRANGEMENTS 


Using the procedure | graphs of momentary angular speed 
of power train parts have been prepared which schematically 
are presented in Fig. 3. Corresponding results are presented 
in Fig. 4. 

As it can be noticed, the quantity that varies most intensively 
is the shaft angular speed of the engine represented by the ©, 
mass inertia moment in Fig. 3. 

The concept of vibration limitation on hybrid arrangements 
consists in assumption that the rotor of electrical machine 
plays the role of a flywheel of the O2 inertia moment. On such 
a system, presented schematically in Fig. 5, it is possible to 
extort a torque impulse from the electrical machine that could 
prevent the variability of its rotor angular speed. 

Modern solutions of electric machine control allow to 
generate electric impulses affecting torque according to a given 
concept within a certain range. A number of ways have been 
investigated how to prevent the changes in rotational speed of 
electric machine rotor. Definitely the most effective way is to 
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Procedure 1. Fragment of PASCAL procedure for calculation of coupling torque moms[j,i] and angular speed om[],i] of drive train elements 


for i:=0 to iii do 

begin 
moms[1,i] 
moms[2,i] 
moms[3,i] 
moms3[i] 


:=(b[4,i] 
:=(b3[i] 


:=(b[2,i]-b[1,i])/e[1]; 
:=(b[3,i]-b[2,i])/el2]; 
-b[3,il])/e[3]; 
-b[3,i])/e3 ; 


om[1,i+1]:=om[1,i]+( + mom[i] + 
+moms[1,i]-tlum*(om[l,i]-om[2,i]) 
)/tet[1]*delt; 


if i>20 then 

begin 
omegasrednie:=0; 
for nn:=1 to 


20 do 


omegasrednie:=omegasrednie+om[1,i-21+nn]/20 


end 


else omegasrednie:=omeg; 


if om[1,i]>omegasrednie 


then deltmo:=-przeciw*momsr 
else deltmo:=+przeciw*momsr; 


moc 


:=deltmo*om[1,i]; 


energia:-delt*moc; 


om[2,i+1]:=om[2,i]+(-moms[1,i]+moms[2,i] 


+deltmo 


+tlum*(om[1,i]-om[2,i]) 
-tlum23*(om[2,i]-om[3,i]))/tet[2]*delt; 


om[3,i+1]:=om[3,i]+( 


-moms[2,i] + moms[3,i] + moms3[i] 
+tlum23* (om[2,i]-om[3,i]) 
-tium34*(om[3,i]-om[4,i]) 

-tlum33* (om[3,i]-om3[i]))/tet[3]*delt; 


om3[i+1] 


:=om3 [i] +(-moms3[i] 


Ttlum33* (om[3,i]-om3[i]))/tet3*delt; 


om[4,i+1]:=om[4,i]+(-moms[3,i] 


— momsr 


+tlum34*(om[3,i]-om[4,i]))/tet[4]*delt; 


b[1,i+1]:=b[1,i]+om[1,i+1]*delt; 
b[2,i+1]:=b[2,i]+om[2,i+1]*delt; 
b[3,i+1]:=b[3,i]+om[3,i+1]*delt; 


b3[i+1] 


:=b3[i] +om3[i+1] *delt; 


b[4,i+1]:=b[4,i]+om[4,i+1]*delt; 


end; 


generate a torque proportional to rotor angular acceleration. So 
the problem is being reduced to the generation by the electric 
machine a torque counteractive to the engine torque. Alas, it is 
impossible by electric motor to generate a counteractive torque 
on the present stage of electric machine control. There are two 
reasons of such situation. Firstly, generation of torque maximal 
over engine operation cycle would require an electric motor of 
mass exceeding many times the mass of combustion engine. 
Seondly, it is impossible to react imediately to the changing 
engine speed. 

Therefore it has been decided to affect the angular speed 
by the electric motor in a simple way consisting in: 
— generation of negative torque equal to the engine average 

torque when rotor speed is higher than average speed, 
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— generation of positive torque of the same absolute value 
over the remaining period. 


The effect of such operation can be evaluated when 
comparing graphs presented in Figs. 4 and 6. 

Disregarding the technical possibilities of so precise 
reaction of electric motor to changing rotor speed, the 
correction of engine speed fluctuations using hybrid systems 
are explicitly limited. The carried out calculation of the real 
level of speed irregularity 5 proves that correction gives over 
double improvement of 6 factor. Above all such operation 
gives the reduction of engine generated noise. Obviously, 
the maximum values of coupling torque in individual parts 
of drive train reduce themselves as well. The values of 


tet1 /tet2 /tet3 /tet4 11/12/13 
0.025/0.113/0.098/0.370 93/10 (20 
tu12/t123 /t.354 /tl3 0/0.0/0.0/0. 


Omega [rad/s] 35R010 tiii227.4[ms] 


20 


D ak nme 
| esq P 


120 160 200  240Time [ms] 


Fig. 4. Change in angular speed of drive train individual parts; 
description in text 


Fig. 5. Schematic of car hybrid drive: I — combustion engine, 2 — fuel tank, 
3 — electric machine operating as motor or generator, 4 — accumulator 
of electric energy, 5 — electric machine and electric energy accumulator 
coupling control, 6 — differential, 7 — wheels 
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Fig. 6. Change in angular speed of drive train individual parts with speed 
correction performed by electric motor operating as a flywheel 


coupling torque on individual sections of drive train have 
been presented in Fig. 7 without correction and in Fig. 8 
— with correction. 
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Fig. 7. Course of coupling torque between individual parts 
of drive train for example presented in Fig. 4 
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Fig. 8. Course of coupling torque between individual parts 
of drive train for example presented in Fig. 6 


The presented considerations about the possibility of 
reduction the torsional vibrations in hybrid systems relate to 
the specific parameters of engine operation. As the additional 
analyses prove, the effect of vibration reduction not always 
is so spectacular as those presented in Figs. 4 and 6 to 8. 
The presented example relates to the case of resonance of 
the highest form (number 3 in this case) and 3 vibrations 
go to one cycle of engine operation. Lower form vibrations 
are dumped less. However, lower form of natural vibrations, 
especially the first form usually happens below the engine 
idle run speed, so it does not cause negative effects within the 
range of engine run. 


SUMMARY 
The problem of torsional vibrations in drives incorporating 


reciprocating engines is still considered as difficult for 
satisfactory solution. In particularly on direct injection diesels 
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what is now a rule, double or even triple inertia flywheels are 
used in order to reduce torsional vibrations and accompanying 
noise. In extreme cases the use of viscous dampers becomes 
indispensable. All these solutions contribute to an increase in 
engine manufacturing costs, increase its weight while the use 
of viscous damper leads to energy dissipation which results 
in elevated fuel consumption. As it turns out by the efforts 
on application of hybrid drives to car drive trains a chance 
emerges to limit vibrations and noise by the use of hybrid 
drives primarily provided for ecologically justified coupled 
combustion and electric drive. The solution proposed in present 
study creates a chance to make use of instantaneous energy 
surplus in two aspects: 


— energy of car braking is accumulated for its use during 
acceleration when further ride is possible 


— energy relative to the increasing shaft momentary speed 
during expansion stroke of combustion engine is 
accumulated for stabilization of rotational speed during 
auxiliary strokes. 
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Problems of welding in shipbuilding - an 
analytic-numerical assessment of the thermal 
cycle in haz with three dimensional heat source 
models in agreement with modelling rules 


Part I 


Theoretical basis of modelling and an analytical assessment 


of heat sources models 


Eugeniusz Ranatowski, Prof 


University of Technology and Life Science, Bydgoszcz 


ABSTRACT 


At the beginning of this paper a short characteristic of the methodology of classic Rosenthal- 
Rykalin solution of temperature distribution during welding is provided. 
In the further part, the requirements concerning process modelling were presented, 
particularly with thermal processes taken into consideration. 
Finally, the Cylindrical-Involution-Normal and Double-Ellipsoidal heat source models 
are presented. 


Keywords: 


INTRODUCTION 


Classic solutions of Fourier-Kirchhoff (F-K) partial 
differential equation, describing heat distribution and leading to 
analytic definition of temperature fields during welding process, 
were first started by works of Rykalin and Rosenthal. The 
temperature fields and other related parameters were found with 
assuming the following simplifications and presumptions: 

a) the analysis is limited to quasi-stationary state, stabilised 
temperature field 

b) heat sources are described by: point, line or plane models 

c) the material physical parameters are constant, temperature 
invariant 

d) welded parts structure is homogenous with isotropic 
physical properties 

e) there are no other inner heat sources during welding 
process 

f) three basic geometrical bodies are introduced in respect to 
the geometry: 

- the semi-infinitely extend solid 

- the infinitely extended plate 

- the infinitely extended rod. 


The assumptions above were some necessary simplifications 
to get the solution but on the other hand affected the preciseness 


of temperature field estimation. The necessary limitations are 
result of methods used for analytic calculation. The adopted 
methods mainly base on integral transform methods and 
Green's functions. 

It should be said clearly that most of heat source analytical 
models as also solid material models are not adequate to 
presently used welding technologies, e.g. plasma stream, laser 
or an electron beam. In the aftermath ofthis, three-dimensional 
heat source model is introduced and also new relation Heat 
Source (HS) — Welded Material (WM) model is created. 

Substantial imperfection of classic Rykalin's and Rosenthal's 
solution is a lack of taking into consideration the heat exchange 
at boundary surfaces. Another important problem is considering 
thermal process as linear, although in fact it's not linear. 


BASIC REQUIREMENTS CONCERNING 
MODELLING RULES 


Modelling is the process of fixing a computational model'. 
The bottom line of any correct modelling procedure is defining 
the parameters characterising the ideology and process 
dimensions. 

An estimation of physic phenomenon during welding 
process, practically results in examining reciprocal relations 
between extensive and intensive parameters. 
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The transport process of extensive magnitudes requires 
estimation of intensive parameters during welding and is 
performed by using such procedures as transient Lagrangian 
or steady state Eulerian formulations of thermal cycle. We 
define an Eulerian (moving ) frame with origin at the centre 
of the source and co-ordinates (x, y, z) — Fig. 1. For Cartesian 
co-ordinate system (X,, Y,, Z,) which remains stationary for all 
time t and loading history, a Lagrangian co-ordinate reference 
is defined. Suppose a heat source is moving at a constant speed 
in the positive x, direction. 

The transformation from (x,, Y» Z,) to (x, y, Z) is given by: 
X=X,-Viy=y,,Z=Z,. 


Zo Z 


Fig. 1. Lagrangian and Eulerian co-ordinate reference 


Specified phenomenon model is connected with presumption 
that the model represents the physic reality well enough. In 
addition, the characteristic feature of the model is ability to replace 
the actual object with the model, still committing transmitting the 
experimental results into actual research object. 

Definition of theoretical structure of research object is 
performed with use of: 

e physic model, describing the actual object 

* mathematical model, being an equation or system of 
equations, describing processes together with boundary 
conditions, characteristic for given phenomenon. 


Above mentioned statements are fundamentally in 
compliance with Fourier’s observation — who stated in “Théorie 
analytique de la chaleur" that “...every phenomenon connected 
with heat depends on very little, general facts the way that every 
problem can be formulated using the language of mathematical 
analysis..." [2]. 

The basic physic laws are described mathematically, 
receiving first generation of mathematical models which are 
last modified (some of those modifications have empirical 
connections). In compliance with Riemann consolidations, 
real basic rules like heat flow, can occur for infinitely small 
magnitudes and must be defined as partial differential equations. 
Integration of such equations leads to settlements and rules for 
variables of space and time e.g. estimation of temperature 
distribution for a plate of optional thickness during welding. 
High conformity and quality of anticipating 1s obtained 1f 
definite phenomenon modelling rules are fulfilled. 

Correct executing of research requires using the proper 
models and conditions with dependence on their likeness. 
Knowing the equations describing given process, for example 
heat flow according F-K equation the criterions and constants 
of similarity can be found. 

If a given variable is marked x, in actual object and x; in 
model object then the connection can be created: 

X X. X 
ie eke = (1) 

Modelling scale is a transforming factor. Therefore the 

following relation can be given: 
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Dt) O(GE L6 0 gig) (2) 
For similar phenomenon there is a relation given: 
DDR OD) 
Equation (3) is fulfilled if c, is in conformity with: 
D(C, 0,040 )=1 (4) 


D(C, Cys C5500, ) = 


Di pM 0s NiO Ris (5) 


DOG X3,X35 Xp) 


In compliance with Federman’s settlements it’s known 
that if a number of variables used in experiment is equal to n 
and a few of them are independent k < n, then there are (n-k) 
equations expressed by modelling scale c,. Therefore a few of 
dimensionless numbers may be proposed by an experimenter and 
be the constants of similarity. The remaining ones are described 
by dimensionless numbers and Federman’s settlements and are 
called similarity criterions. The indispensable conditions for 
similarity existing between the models and actual objects is 
describing the happening physic processes by: 

e the same differential equations with appropriate boundary 
conditions 
e the similarity criterions verification. 


Let’s take into consideration F-K partial differential heat 
flow equation as a starting point — for stationary co-ordinates 
system X, Yo; Zo: 


. oT 
div(à -grad T)-c, - por = —q, (X0; Yo Zot) © 


The classic Rykalin’s and Rosenthal’s solutions assume that 
physic parameters describing the process model are constant 
whereas they are in fact non-linear for real welding system as 
a result of dependence on temperature. 

According to that we have the following characteristic 
[following equations (1) — (4)]: 


Process model Actual welding system 


À = const MT) # const 
c, = const c (T) # const 
p = const p(T) # const 
ACT c (T) T 
E o ee P. 3 (7) 
À C, p 
c Z6, £706, 
D(c,, c, c) £1 (8) 


Moreover, q, doesn't fulfil geometric similarity condition 
for majority of contemporary welding methods like arc welding, 
GTA, GMA, plasma stream and others, assuming point, line or 
plain heat source model 

Considering possibility of universal shape modelling and 
to fulfil conditions (4), (5), two main models are considered: 
Cylindrical-Involution- Normal (C-I-N) [3] and The Double 
Ellipsoidal Configuration of Source (D-E) [4]. Both models 
represent effective extended abilities of 3-dimensional heat 
distribution. The model that previously was used is the Gaussian 
surface flux distribution model: 


q= SE epeka +y- © 


Unfortunately this one does not reflect heat power input 
change for “z” variable. For example Wei and Shiar? constituted 
that under high-intensity laser beam welding the incident 
energy rate distribution is assumed to be Gaussian distribution 
and HS model is idealised by a paraboloid of revolution. This 
is one of main reasons why two more accurate models are 


implemented. 


THE CHARACTERISTIC OF HEAT 
SOURCE MODELS 


The mathematical expression of the cylindrical-involution- 
normal C-I-N model is: 


q,7 qux: Exp k(x’ + y’)-K,-z]-[I- u(z—s)] (10) 


66,99 


As it is seen “z” variable is not squared in contrast to other 
variables. This is an intentional purpose of this model to let it be 
easily integrated when using integral transformation method. 

D-E model doesn’t have this feature which makes it difficult 
to achieve precise solution while transforming. 

Let's calculate q qax factor for C-I-N. Provided that the 


following condition must be preserved: 
o0 o0 o0 


Q= | J fasmax exp E k(x? +y*)-Kz-z]- 


0 —oo— oo 
:[1- u(z- s)]dxdydz 


and after integrate operation upon several variables we 


obtain: 
oo oo Ja 
Q= | | Fe dec xP C ky KZ) iy 


0 -% 


(11a) 


-(1- u(z-s)dydz 


dmax EXPE K,:z)[l-u(z-s)]dz (11c) 


T -7 
= ——— qymax’ EXPO K, -S)———— qymax(11d) 


k-K, k-K, 
thus: 
Q 
= —___+___—_.-k -K 11 
IST NI 
and finally: 
Q 
ge——— — kK, 
n -[1- exp K , -s)] T 


:exp[- k(x? t y?) -Kz:z] [1 -u(z -s)] 


The [1-(u(z — s))] factor illustrates that q, = 0 if z >s, see 
Fig. 2. 

By changing: s, k and K, factors, C-I-N model can represent 
all presently used heat sources. 

Let's use equation (12) to analyse the shape of the surface 
that limits the volume of effective C-I-N affect. Surfaces of 
constant values can be obtained by comparing: 


y= quax exp[-k(Gc y^) -Kz:z) 
-(1- u(z-s))]= const 


(13) 


1 - u(z-s) = u(s-z) 


S Z 
Fig. 2. Heaviside s functions 


where: 
k-K,-Q 


vax 1-[I- exp(-K, -s)] 


(according to (11e). 
Presuming z < s we receive: 


exp[-k - (x? + y?) -K,z]- SOME A 
QvMAX 


-k:(x2« y5)- K,:z-1n(A) - -B 


a B-k-(x°+ y) 
K 


Solution (14) gives a family of paraboloids of revolution. 
Taking into consideration [1-u(z-s)] factor, it is obvious that 
some of the paraboloids may be “cut”. A few examples of 
paraboloids are shown in Fig. 3. All of them have the same 
height h = s. 


(14) 


zZ 


AETHER * 
a IIT NN 


Y 


Fig. 3. Examples of the paraboloids of constant values of q, 
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The paraboloid that limits the volume of effective input 
affect can be obtained using Rykalin’s suggestion that HS 
power loses its effectivity when: 


2 2 
dymaxexp[-k- (x^ y) -K,: z] = 0.05q max 
-—k-(x°+y9)-K,-z=3 
3 k 
geo 2uafy? xy?) 
K, K, 

Equation (15) is an expression that describes effective input 
affect volume. In dependence on h = 3/K, value, the mentioned 
paraboloids can be more or less slim — see Fig. 4. At extreme 
conditions one can obtain several, well known heat source 
models like: 

a) for K, — 0, k # 0 — cylinder with height = s — disc heat 
source 
b) for K — oo, s— 0, k # 0 — surface heat source. 


Z£S (15) 


b) 
I5 c 
A r 
rı ro = 0 
B 
roS r< T2 
z A B C 


Fig. 4. The paraboloid limiting C-I-N effective input power volume (a) 
and power in dependence on distance from centre and depth 


Possibility of changing K, ,k, s and Heaviside's function 
u(z-s) values makes C-I-N heat source a universal one. 

The Double Ellipsoidal heat source model (D-E) (Fig. 5) 
established by J. Goldak [4] et al. has the following form: 


q D £630 . 
$ a-b-c-n- yn 


_ £64520 (16) 


Tam A i 
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Fig. 5. The Double Ellipsoidal configuration of Heat Source model 


D-E arrangement of volumetric heat sources is used for 
moving welding heat source of deep — penetrating surfacing 
or butt welds. The front half of the source is the quadrant of 
one ellipsoidal source, and the rear half is the quadrant of 
another ellipsoid. According to Godak [4] et al. some peculiar 
big pressure should be exerted on models, which are the 
combination of two or more usual heat source models, in order 
to effectively reflect the actual welding process. The structure 
of equation (16), which is an excellent one for numerical 
estimation of thermal cycles, has a few inconveniences when 
used for about finding analytic solution from methods described 
in this investigation. However, some estimation can be done 
as it is worthwhile to make some efforts, in order to obtain the 
necessary solution. 

The adopted C-I-N model is favourable for simulation of 
welding process of high concentrated energy source used, like 
laser beam welding or electron beam welding. This one may 
also be used for the simulation of arc welding process. 


CONCLUSIONS 


Base on the fact that substantial imperfection of Rosenthal's 

and Rykalin's solutions relies on: 

* considering thermal process as linear 

e lack of taking into study the heat exchange at boundary 
surface, 

most of heat source analytical models and solid material models 

are not adequate to presently used welding methods, there are 

established the basic requirements which concern modelling 

rules of welding process. 

It is constituted that the indispensable conditions for 
similarity existing between the actual object and model describe 
the happening physic process by: 

* the same differential equations with appropriate boundary 
conditions, 
e the similarity criterions verification. 


In order to fulfill above conditions, at first two heat source 
models are considered: C-I-N and D-E which represent effective 
extended abilities of 3 - dimentional heat distribution. 


NOMENCLATURE 


- temperature, [?C or K] 

- thermal conductivity, [W cm" K"'] 

- specific heat, [J kg! K`] 

- mass density, [kg cm?] 

- time, [sec] 

power input in volume, [W cm?] 

- net power received by the weldment, [W] 

- a factor designating the HS concentration, [cm?] 

- HS penetration depth, [cm] 

u(z-s) - Heaviside's function 

a., b., c,, a, b, c, - ellipsoids semiaxes, [cm] 

- fraction of heat deposits in front and rear 
quadrants 

- power density distribution inside the front 

and rear quadrants respectively, [W cm?] 


2 wo wer 
I 


©: - modelling scale 


i 


K, - involution factor of HS, [cm] 

HAZ - Heat affected Zone 

C-LN - Cylindrical-Involution-Normal heat source 

model 

D-E - Double-Ellipsoidal heat source model 

HS - Heat Source 
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Dynamical control of dimensional quality 
of large steel structures in production 
enterprises of low - level technological support 


Radostaw Rutkowski, Ph. D. 
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ABSTRACT 


In this paper problems of dimensional control in the building process of large steel structures 
are discussed. Results are presented of the research aimed at the improving of dimensional 
quality of prefabricated ship hull structures built in enterprises of medium — or low- level 
technological support. Problems of dimensional control in building process of large steel 
structures are highlighted. A solution is proposed of dynamical control of measurement 
operations as well as their optimization with taking into account current production process 
parameters. An example of calculations carried out in the frame of the presented solution 
is also attached. 


Keywords: shipbuilding metrology, technology, large steel structures, shipbuilding, 
dimensional quality, mathematical modelling. 


INTRODUCTION 


As dimensional quality problems are not very much spread 
in shipbuilding circles, basic introductory information on the 
issue is presented below. 

The main task of a shipyard is obviously to produce ships 
of the highest quality at the possibly lowest cost, i.e. arising 
from economical consumption of materials and the possible 
shortest production cycle. Moreover reaching the economic 
success depends on keeping the assumed contractual parameters 
such as a.o. ship’s length, depth, breadth, minimum keel 
plate deformations, and finally- ship’s speed, deadweight 
and fuel consumption etc. The reasons lead to very high 
quality requirements in the sense of design of technological 
processes and rationalization of processes of manufacturing 
and production quality control. Ship hull is assembled of many 
structural elements. Nowadays built ships are often consisted 
of more than a hundred sections. In Fig. 1 an example sectional 
division is presented. In order to make correct matching and 
joining the hull elements possible it is necessary to keep their 
dimensions, shapes and mutual position in compliance with 
design of the ship. In realizing the aim measurement control 
operations serve as main tools applied at every building phase 
beginning from lofting through cutting, bending, up to final 
assembling the ship hull after launching. They make it possible 
to check compliance of design of a given structure with its 
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realization during production cycle up to its end. The problems 
are covered by shipbuilding metrology. 

Dimensional quality problems of small-size (up to 500 
mm) and medium-size objects (up to 1150 mm), which 
concern both technological processes and means of production 
are fully mastered. For manufacturing such objects the 
GPS system (Geometry Product Specification) containing 
methods of dimensioning, dimensional tolerancing and fitting 
collaborating elements has been developed. Its description 
is given in [2]. In many cases the solutions proposed in 
the frame of the GPS system do not concern large-size 
objects which are produced mainly in shipbuilding industry. 
Scale of the dimensional quality problem is large in production 
enterprises of a technological level similar to that of Polish 
shipyards. As results from the author’s experience in such 
shipyards the increase of labour consumption resulting from 
repair of dimensional quality errors reaches as much as 
30 + 40% of total labour consumption for building ship hull. 

The problems in question are not widely discussed. 
Moreover the published solutions are mainly dedicated to 
shipyards of a high-level technological support. Such level has 
been reached by realization of complex restructuring projects for 
shipbuilding industry at the expense of very high initial outlays. 
To use the then-presented solutions in production enterprises of 
characteristics close to Polish shipyards is very difficult, if not 
impossible at all. Their features are discussed below. 


Fig. 1. Hull sectional division of a ship built by Szczecin Shipyard 


ASSUMPTIONS AND LIMITATIONS 


As already mentioned the dimensional quality problems in 
production enterprises of characteristics similar to shipyards 
are very individual. It arises first of all from the following 
features: 

e high production intensity 

* large changeability of manufacturing process, that makes 
permanent control of dimensions during the process 
necessary 

* seperate dimensional quality requirements for particular 
units 

* low investment possibilities 

e alackofroofing over hull assembly sites (which produces 
significant problems resulting from deformations due to 
insolation) 

* alarge degree of co-operation (many structural elements,e. 

g. hull sections and blocks, are produced by subcontractors 

in various manufacturing conditions). 


The conditions make application of versatile measurment 
systems very difficult and even impossible in most cases. 
Moreover application of electronic measurement instrumentation 
which requires specialty measuring stands as well as significant 
financial outlays is rather troublesome. 

The measuring methods most used in the above mentioned 
cases are those based on optical measuring instruments and, as 
far as conditions allow, electronic tachymeters with application 
of the pole measurment method. The crucial assumption 
(requirement) put for elaboration of solutions presented in this 
paper was to use the measurement methods which have been 
applied so far. 


Moreover the necessity of application of optical measuring 
instruments results from the high degree of co-operation in 
the area of prefabrication of hull sections. In most cases the 
subcontarctors do not have at their disposal production sites 
suitably prepared to carrying out dimensional quality control 
by using advanced measuring techniques. 


IDEA OF DYNAMICAL DIMENSIONAL 
CONTROL OF LARGE STEEL 
STRUCTURES IN PRODUCTION 
ENTERPRISES OF MEDIUM- AND 
LOW- LEVEL TECHNOLOGICAL SUPPORT 


In structurebuilding process measuring processes intertwine 
with technological ones. We deal here with a complex technical 
system, to use system approach is therefore necessary for 
analysis of measurement control operations. Geometrical 
control system in the area of shipbuilding metrology is defined 
by this author as a set of elements of metrological operations, 
their manners and set of couplings between the elements and 
those of the technological process itself. The models in question 
make dynamical control of measuring processes possible on the 
basis of necessary accuracy standards as well as information 
gained from production process under way. The accuracy 
standards are determined by required dimensional accuracy 
of prefabricated and assembled hull structure elements. The 
problem is discussed in the further part concerning verification 
of the results obtained in modelling process. 

The discussed models represent control and measurement 
systems mainly by using mathematical description as well as 
- toa much smaller extent - linguistic one. This is a simplified 
description of the process features considered important with 
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regard to application of the model. In the discussed case the 
task ofthe model is to indicate optimum parameters of control 
and measurement processes (i.e. their system). 
The parameters are first of all the following: 
* measurement accuracy 
e measurement method to be used 
* kind and model of measuring instrument 
* co-ordinates of measuring instrument stand 
* co-ordinates of base points. 


Number of the parameters and their characteristics depend 
on a type of measuring task. 

The overall schematic diagram of the discussed models is 
presented in Fig. 2. 


INPUT MODEL INPUT 
VALUES VALUES 
CONSTANT EQUATIONS VARIABLE 


CALCULATION 


OUTPUT 
VALUES 


OPTIMUM 
PARAMETERS OF 
THE SYSTEM 


NO 


Fig. 2. Schematic diagram of modelling process 


The modelling process begins from detail description of 
a measurement task for which control & measurement system 
has to be designed. Next, the following input quantities should 
be determined: 
* constants characteristic for a given shipyard, available in 
the designing phase 
* variables which charecterize changes occurring in 
production process and affecting a modeled system. 


The successive phases of modelling process are as follows: 
* calculation processes performed on the basis of model's 
equations 
* verification process. 


„The heart" of the system is composed of the model’s 
equations which determine input quantities on the basis of 
which are determined parameters subjected to verification in 
the successive step. The equations are elaborated on the basis 
of mathematical description of measurement methods. In order 
to present the described system it is necessary to highlight 
measurement methods used in shipbuilding industry, and their 
mathematical models. 


MEASUREMENT METHODS 


As already mentioned the control & measurement 
operations during building process of ship hull are usually 
performed with the use ofthe pole method based on electronic 
tachymeters, as well as classical methods (the optical levelling 
method and reference line method) by using optical measuring 
instruments. 

The pole method based on tachymeters is widely known and 
used during building process of large structures. Its application 
to shipbuilding metrology was described a. o. in [4]. And, as 
results from the author's experience, the shipbuilding metrology 
based on the classical optical instruments should be highlighted 
especially as regards its mathematical description on which 
further analyses are based. 
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The main characteristic feature ofthe discussed measurement 
methods is the necessity of carrying out measurements within 
inclined co-ordinate frames as well as in dynamically varying 
systems as ship hull assembling sites are usually inclined 
relative to levelling plane. Moreover a part of ship's outfitting 
operations is carried out on ships lying along berth (floating 
ships undergo oscillating motions). The conditions make direct 
application of typical measuring methods based on constant 
levelling plane determined by means of levels or self-levelling 
mechanisms, difficult and even impossible in many cases. The 
methods in question allow for making the reference frame 
independent of the levelling plane. Below is presented a short 
description of the mentioned methods as well as a procedure 
for elaborating their mathematical models. 


REFERENCE LINE METHOD 


The reference line method finds broad application both 
to measurement of deformations and assembling control of 
objects. It is known in a few versions which differ to each 
other by measuring technique and kind of quantities to be 
determined. Measurements performed by using the method 
mostly relate to two mutually perpendicular planes: vertical 
one containing reference straight line, and horizontal one. 
During control & measurement operations carried out in ship 
building process the two -plane system must be replaced - for 
the above discussed reasons (to make the reference system 
independent of the levelling plane) - with the system of two 
mutually perpendicular planes II, , IL, the first of which is 
parallel to ship's base plane (Fig. 3). 


Fig. 3. Reference line method 


The reference straight line is usually formed by sight line 
of theodolite or levelling instrument. The method in question 
is used to determine the coordinates y, z,; the coordinate x, 
plays an auxiliary role (as a calculation parameter). For every 
quantity other measuring technique is used. It results from 
that the quantity y, is related to the plane IL, whereas z, to II. 
Hence seperate mathematical descriptions were elaborated 
for determining the coordinates y, and z, The coordinates z, 
describe values of the sections perpendicular to the reference 
line and located on the plane IL. Fig. 4 was used to elaborate 
the mathematical model for determining the quantity Zoo 

In the points E, F, and P measuring rods are placed 
perpendicularly to the II. The sight line of measuring 
instrument determines readouts ofh,, h, and h, The searched- 
for value of z, is determined by using the following relation: 


Ze = PP" Lg (1) 
On transformation the following is obtained: 
hp-h 
Xp 


ing instrument 
ght line of measuring, jnstrur 


hg 


Fig. 4. Reference line method 


The above given formula constitutes the mathematical 
model for determining the quantity z, by means of the reference 
line method. 

To carry out calculations by using Eq.(2) on the assembling 
stand is labour-consuming. To shorten them the measurement 
executor tries to reach the situation when h, = h, In such case 
the determining of the quantity z, is simplified to the following 
form: 

2,-h,—h, (3) 
that improves the measuring process considerably. 

The coordinates y, describe the sections perpendicular to 
the reference line and parallel to the plane II,. To determine the 
quantities, an instrument whose sight line is able to build the 
plane II, perpendicular to II, and directed along the reference 
line EF, is necessary. The searched- for value is equal to the 
distance of the point P from the plane II, (Fig. 5). 


Fig. 5. Reference line method 


Fig. 5 was used to elaborate the mathematical model for 
determining the quantity y, The model in question was formed 
on the basis of the expression describing the distance of the 
point P(x, lp Zp) from the plane IL, which contains the points 


EX p Yp Z,), F(x, Yp Z,), G(x, Vo Zo): 
Rip +Sx, +Tz, +U 


Yp (4) 
JR +9 +T? 
where: 
R=(x,°Z, XZ X, ZORXSUASCHXOUAS Ke Z) 


S=( Ye 2p Yp’ Zo Yo ` Ze t Yg ` Zo Y Yp ' 4 T Ya’ Zp) 

Teu Xo Yr’ Xg t Ya’ Xg T Yg’ Xp Xp’ Ya Yr Xo) 

U=Cy, "Xe Ze Ye Xp Zo T Ya’ Xp’ Zg t Yg Xe zt 
-Yp ' Xp UN ey "X Zp) 

RI, + Sx, + Tz +U >0 


On the assumption that the auxiliary coordinate system 
Oxyz associated with the planes II, and II, satisfies the 
following conditions: 

e x — axis of the system forms edge of intersection of the 

planes II, and II, 

* y-axis lays on the plane II, 
e z- axis lays on the plane II, 


the following can be assumed: 

* WYT YO 

° Zz-z,-0 (5) 
e X D 


On insertion of Eq. (5) to Eq. (4) the following was 
obtained: 


ys 74 (6) 


As results from Eq. (6), in the assumed system the 
measuring of the quantity y, is equivalent to the reading-out 
of value of the quantity 1, on the measuring rod put against the 
object in the point P and placed perpendicularly to the plane IL. 
The compexity degree of Eq. (4) makes calculating y, - values 
greatly difficult. In practice one tries to reach the situation where 
the conditions (5) are satisfied that makes using Eq. (6) for the 
calculations possible. 


OPTICAL LEVELLING METHOD 


As already mentioned in shipbuilding metrology it is 
necessary to relate measurements to planes or near-planar 
surfaces inclined relative to levelling plane by the small angle 
©. In such cases appropriate measurements are performed with 
the use of the optical levelling method in which the reference 
level is not used, but the auxiliary plane II, is applied (Fig. 6) 
relative to which measurements are carried out. The method 
may be applied both to carrying out measurements on land- 
based assembling stands and on floating ships (floating dock, 
ship after launching). It is schematically presented in Fig. 6. 


4 


GIULI. a A 


UE, 


A 


Fig. 6. Optical levelling method 


The crucial task in conducting the measurements by means 
of the optical levelling method is to determine the coordinates 
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Z, of the points P, in an orto-cartesian coordinate system Oxyz. In the case when the plane II, constitutes the reference surface, 
the axes x, y of the system are assumed to lay just on the plane, giving them this way orientation appropriate to given purposes. 
The plane II, is determined by three suitably chosen points A, B, C which usually lay on a measured structure. The points are 
called the base. 

Value of z, is calculated in accordance with the following relation: 


zd (7) 
where: 


Z, —.Z'- coordinate of the point P, 
d . 


Zp © Zp 


Value of z, n be determined by using the relation which expresses the equation of the plane II, which contains the points 
AGC, Vo h,), BZ, yp h DE C(x, Vo h,): 


Xp Yp Zp, I 
x h 1 
A Ya A -0 (8) 
Xg Yg hg ] 
Xo yc he 1 
On transformation, Eq. (17) can be presented as follows: 
-— (AX syp + AXgpyA— AX,pyp)hc + (AX ap¥c — AXAcyp- AX copy a p+ (Axgcyp e AXgpyc + AxXepyg Ny (9) 
Po 
AX Ve PAR aca Pag 
where: 
AX yp 7 X4 XX 
AX yp =X,—X, 
AX,, =X,—X, 
Ax,c = Rey 
AXc, CX Xe 
Axe = Xe T Xg 
On insertion of the right-hand side of Eq. (9) to Eq. (7) the following was obtained: 
de (Ba FAXgy,— AX ap Vp het (NX ep Vo AK ac ¥p— Axes Np t Axa Axg yo Axon Vp Ba h 
P— P 
AX Vo FAX5oYa AXaoYa (10) 


This 1s the mathematical model of the optical levelling 
method. 

The carrying out of calculations on the basis of the above 
given expression makes the control & measurement process 
much longer.The calculations in question become simpler when 
h,=h, =h,. In such situation Eq. (10) takes the following 


form: 
z,-h,-h, (11) 


In practice one should aim at the situation when to 
apply Eq. (11) is possible as it greatly simplifies calculation 
process. 

By analysing the above presented mathematical descriptions 
of the measurement methods it is easy to come to the conclusion 
that accuracy of determining measured values is first of all 
affected by errors of readouts taken from measuring rods. The 
quantities are described in the subject-matter literature, e.g. [3], 
by means of some approximate formulas which are decisively 
too-low exact for purposes of shipbuilding metrology. Therefore 
the quantities were determined experimentally. To this end 
a special testing device was elaborated by using which tests were 
performed for selected measuring instruments. In consequence 
of the tests regression equations which describe values of the 
errors in question, were obtained depending on aimed length 
values. The equations were elaborated individually both for 
particular instruments and measuring stands. The obtained 
results confirmed the above given thesis on too-low exactness of 
the approximate formulas describing the errors in question. 
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A way of making use of the obtained equations has been 
presented with the use of a calculation example given below. 


MODEL EQUATIONS 


During modelling process calculations are carried out on 
the basis of model equations elaborated by using mathematical 
description. The models were used for elaborating accuracy 
description of measurements performed by using particular 
methods. The description was prepared with the use of the 
principle of proliferation of mean errors. As a result of the 
operation the following was obtained: 

For optical levelling: 


2 2 2 
N Q 
PONES a 
BUT ah, oh, 


2 2 
Oz " Oz. i OZ 
mi, +H — | mx, +| — 


OX. 


da: 2 . Oz, 2 : n. 2 : (12) 
+ — | mi mjs +| —— | mk. + 
OX, OX 


n — number of points P, 
m,,m,,m,,m, — standard deviations of respective 
€ Pi readouts from measuring rods 

m, m d. 

m4,m “4 m 8 
Yg! Xe Ie sati inati 

m, > m“, m, — standard deviations of determination of 
c P P: 


coordinates of the points A, B, C and P.. 


To the variable z, the parameter ,,i” was introduced, (z, ), 
which describes a number of the points P whose coordinatès 
»Z are to be determined in one measuring process. 


For the reference line method the following was 
achieved: 


Forz: 
E 2 2 
Z. 
= "a nt —— ms, + 
E B (13) 
Oz, Y Oz, | Oz, Y 
Z Z Z 
+ | my +) — | m’x, +} —* m. 
ch, OX, 2 
where 
X, X, — coordinates of the points E and P 
h,, hp h, —  readouts from measuring rods placed in the 
points E, F, and P 
mm, — standard deviations of determination of 
d coordinates of the points E, P 
m,m, >m, standard deviations of appropriate readouts 
i from measuring rods. 
For y, 
y yy | o ,( Oye) 
m$, E jm 3 a. == 
Kr Oyg Zg 
oy 2 Di 2 oj 2 
+ iJ ins, + p m$, + - mz,+ 
Z 
2 ay 2 ay 2 
f - mis +| —— | my, + —— | mz 
Ka Ya Za 
Bye) > ,( Qe) 2 ,( Se) 
+)? | mx, +] =? | mi, +| =* | mz, 
ox, al, Oz, 
where 


— standard deviations of determination 
of appropriate quantities. 


Mathematical model of the pole method as well as 
description of accuracy of measurements carried out by using 
the method, was presented in [4]. 


CALCULATIONS AND THEIR 
VERIFICATION 


Description of procedure of using particular equations 
would require a very extensive elaboration. Therefore the author 
considered it reasonable to present only a selected calculation 
process based on the equations of the optical levelling method. 
Calculation procedures for the remaining equations are carried 
out by using similar principles. 


The calculations consist in finding optimum 
parameters (decision variables) for assumed input data. 
To this and, appropriate optimization methods were selected 
to take into account complexity of particular models. In the 
case in question searching for optimum parameters consists 
in finding an extremum of the function (15) which describes 
standard deviation of a given measuring process. In the case 
of measurements carried out by means of the optical levelling 
method the function takes the following form: 


F2». — min (15) 
il 
171,2,..n 
where: 
n -— number of the points P (measured ones), 
r, — quantity which expresses determination accuracy of the 


Pi coordinate ,,z” for a given point P. 


The task is realized by finding a value of the vector u, 
Eq.(16), which minimizes the above given function. 


u= By Y ap Yp Xc, Yc, Xp Yp] (16) 
where: 
xi Ya Xp Veo Xo Vo — coordinates of the base points, 
X,» y, — coordinates of the measuring stand. 


The quantity r, was defined on the basis of Eq. (12). Its 
form is given by thé following expression: 


2 2 2 


Oz OZ OZ 
li. = Biel, ^ + Pi «ti + Pi «ti + 
: | Oh, oh, oh, 
2 2 2 
z OZ OZ 
+ Pi t; Pi + Pi + 
Oh, * Y OX. Oy 4 
: (17) 
OZ, OZ, OZ, 
+ -| +| —— |I + i | + 
i; Oy, Oxo 
2 
+ dz, dz, dz, 
Oyc OX, OY p, 
171,2,..n 
where: 
n — number of the points P 
t,, tẹ t, t, — equations which describe standard deviation of 


i 


readouts taken from measuring rod. 


In the above given equation, values of standard deviations 
of determination of coordinates of particular points (see Eq. 12) 
were not taken into account as they can be considered constant 
and mutually identical. The performed calculation tests which 
consisted in inserting the constants to the optimization equation, 
showed that they did not affect the obtained results. 

On the basis of the analysis of values of readout errors the 
following general form of the equations ,,t” which desctribe 
standard deviation of readouts taken from measuring rod, was 
assumed: 


t, - ag; ^ bg, +c CE) 
where: 
1=A,B, C,P, — location points of measuring rods, 
a, b, c — coefficients of the equation which describes 


aiming error, 
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sight line length for readouts 
taken from measuring rods 
located in the points A, B, C, P, 
coordinates of the points A, B, C, 
P, as well as the measuring stand 
location point J, respectively. 


gi -Je-x) (y; - 


Xp X» y» y; ~ 


Values of the coefficients a, b, c of the equations ,,t” result 
from a regression analysis performed for the obtained test 
results. Example results obtained for measurements carried 
out on slipway are presented in Fig. 7. 

Ni-021A levelling instrument 
y =0,0001714 - x^ + 0.0017857 - x + 0.3242857 


12 
11 
10 
09 


aiming terror [mm] 


0 10 20 30 40 50 60 70 80 
distance [m] 


NI 004 levelling instrument 
y = 0.0000548 - x? + 0.0039048 - x + 0.2385714 


aiming terror [mm] 


"0 10 20 30 40 50 60 70 80 
distance [m] 


Fig. 7. Aiming error and reading-out error 
of values taken from measuring rod 


Values of the searched-for vector u, Eq. (16), are to satisfy 
a given set of equations or inequalities — constraints. In the 
discussed optimization process the constraints result from an 
allowable location of measuring stand and base points. The 
location can be determined as an area of set of coordinates. 
If this is an area, the inequality constraints a < c < b are used, 
where c stands for the coordinate x or y of the points A, B, C, J, 
and a, b consitute its lower and upper constraint. And, if the set 
of coordinates is selected, the equality constraints: c= (d, ..., d ), 
where d, stand for values of the above mentioned coordinates, 
possible to be selected, are applied. In practice, location of base 
points is usually described by pairs of coordinates, but location 
of measuring instrument is indicated as an area. 

In shipyard practice the coordinates ,,z "are determined even 
for a dozen or so points P. It generates a dozen or so objective 
function components described by Eq. (17). Finally for the 
case in question a very complex formula has been obtained. 
To make applying the proposed solutions to shipbuilding 
metrology possible it is necessary to use a software by which it 
would be possible to fast elaborate and minimize the discussed 
functions. 

Choice of such software and optimization method is rather 
subjective, depending on which computer and software is at 
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hand. This author selected the Mathcad software. The choice 
was supported by multiple test calculations. The software offers 
two optimization methods: 

e quasi-Newtonian one 

* method of coupled gradients. 


Calculations were performed by using both the methods. 
Differences in their results were found small. Better results (the 
objective function reached smaller values) were achieved by 
using the quasi-Newtonian method. Moreover superiority ofthe 
method was supported by that some problems with satisfying 
constraints arose when the method of coupled gradients 
was used. Summing up, the Mathcad software satisfies the 
following basic requirements for solving the calculation tasks 
of the kind: 

* repeatability of results 

e fulfilment of constraints 

* short time of calculations (if suitable calculation modules 
are prepared in advance the time of realization of a complex 
optimization task by using the optical levelling method will 
not exceed one hour). 


The detail description of the software as well as its 
optimization methods can be found a.o. in [1]. 

In order to better highlight the described process in the 
further part of the paper a simple example of the modelling 
process which covers also the optimization process based on the 
optical levelling method, is given. Of course there is possible 
to present more complex calculation processes. However it 
would require to make a seperate elaboration. 

The next step of the modelling of a given system is the process 
of verification aimed at selection of parameters ofthe system, which 
make it possible to reach such accuracy of measuring process that 
satisfies technical and economical requirements. The technical 
requirements are associated with manufacturing tolerances on the 
basis of which minimum accuracy of measurements is determined. 
However the economical requirements are connected with cost 
of performance of measurements as well as accompanying 
operations (e.g. stopping production process, reorganization of 
a production stand on which measurements have to be carried 
out etc). The latter criterion is deemed secondary for the process 
in question as that orientated towards fulfilment of technical 
requirements is of priority importance. 

In order to determine minimum accuracy to apply the 
commonly used ,, golden principle of measuring technique” of 
Berndt, consisting in that standard deviation of a measurement 
should not exceed 10 % of manufacturing tolerance, was 
proposed. And, in the case of obtaining a measurement result 
placed within uncertainty area (see Fig. 8) a technological 
analysis of production process under way should be performed. 
In the situation when the procedure indicates necessity of 
increasing measurement accuracy the modelling should be 
repeated for a modified technical criterion and to perform again 
measuring process. 

In the case when obtained results do not satisfy verification 
criteria, input constant data should be changed by applying 
another: 

* measuring instrument 
* measuring method. 


As already mentioned, the calculation and verification 
processes make it possible to determine the optimum parameters 
of modelled measurement process, such as: 

* measurement accuracy 

* measurement method to be used 

* type and model of measuring instrument 

* coordinates of measuring instrument stand 
* coordinates of base points. 


non-compliance area 


value of dimension 


uncertainty band 
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Fig. 8. Influence of measurement uncertainty on assessment of obtained 
results, where: m e standard deviation of determination of R — dimension; 
T — manufacturing tolerance, compliance area — tolerance area lessened 
by standard deviation of determination of R — dimension; non-compliance 
area — areas outside the tolerance area increased by standard deviation 
of measurement; uncertainty band — the interval within which neither 
compliance nor non-compliance can be stated 


EXAMPLE MODELLING PROCESS 


To better highlight the described processes and their 
possible applications an example of modelling process based 
on a real measurement task is presented below. 

Inthe given example the modelling equation was subjected 
also to maximization process that showed a scale of possible 
profits available due to optimizing process. 

The four phases can be distinguished in the modelling 
process in question: 

1* phase: Determination of input parameters (measurement 
method, measuring instrument, available coordinates of 
base points, admissible location of measuring instrument, 
manufacturing tolerance of an object to be measured); 

2™ phase: Determination of optimum location of measuring 
instrument and base points [acc.Eqs (15) (17) (18)]; 

34 phase: Calculation of standard deviation of ,,Z" — coordinates 
of measurement points [acc. Eq.(12)]; 

4" phase: Verification. 


The designed control & measurement system is aimed at 
checking waviness of inner bottom of a bottom section of ship 
hull during its assembling on slipway. The schematic situation 
diagram of the object to be measured is presented in Fig. 9. 
Admissible area for setting up 
the measuring instrument 


— Admissible locations for base points 


P — Measurement points 


Fig. 9. Schematic situation diagram of the discussed measurement task 


As assumed, the measurements were performed with the use 
of Ni-021A levelling instrument (a relevant value of readout 
error taken from measuring rod is given in Fig. 7). Tolerances 


for measured quantities, established by ship design office were 
equal to + 7 mm. 
In order to keep this presentation as clear as possible only 
final results of particular calculations are attached. 
In the case in question the objective function is composed 
of six components r, : 
Pi 


Fer EL SS a, FI E 
P 


2p, Zp 2p, Zp; 2p, 


(19) 


Particular elements of the objective function are described 
by using Eq. (17). The component t described by means of 
Eq. (18) took the following form for the investigation results 
shown in Fig. 7: 


t, = 0.0001714[«, —x,) +(y,-y,)]4 


+ 0.0017857 J (x, — x) + (y, - y;) + 0.3242857 


The obtained optimization results, i.e. the coordinates of 
the base points as well as measuring stand, are presented in 
Fig. 10. Values of the coordinates are given in meters. In Fig. 
10 the most unfavourable variant is also given. 

a) VARIANT 1 (optimal) 
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Fig. 10. The optimum variant of the measurements, (a), 
and the most unfavourable one, (b) 


For the parameters presented in the figure was determined 
standard deviation of ,,Z" — coordinate of the points P , P,, P}, 
P,, P, P.. For the calculations Eqs. (12) were used. The value of 
the standard deviation of determination of location coordinates 
of measuring rods was assumed to be 5 mm. The obtained 
results are graphically presented in Fig. 11. 

Taking into account the value of manufacturing tolerances 
of the measured object, equal to + 7 mm, one can state that 
the obtained results of measurement uncertainty satisfy the 
basic requirements according to the Berndt principle. When 
after completing the measurements the obtained results are 
placed within the uncetainty area (see Fig. 8) the measurements 
should be repeated with changed input parameters. In the 
case in question the most suitable operation would be to use 
a measuring instrument of a smaller error of readouts taken 
from measuring rod. 

On the basis of the obtained results it can be stated that this 
is choice of location of base points and measuring instrument 
stand which is crucial for accuracy of measurements performed 
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Bl Variant 1 


Standard deviations [mm] 


Pl P2 P3 P4 P5 P6 
Fig. 11. Values of standard deviation of determination 
ofz” — coordinates of the points P, 


with the use of the optical levelling method. The results 
obtained from the above presented calculation process for the 
most unfavourable variant as well as the optimum one reveal 
positive consequences of application of optimization methods. 
The presented case is rather not complex, for more complex 
measuring processes profits due to application of the discussed 
methods are greater. 

Asthe presented example confirms, due to the application of 
modelling process effective control of manufacturing tolerances 
can be ensured by using optical measuring instruments. 


CONCLUSIONS 


* The solutions presented in this paper result from 
investigations in the area of shipbuilding metrology, carried 
out for the sake of production enterprises of medium- or 
low- level technological support. The fact has determined 
the guidelines as well as limitations which appeared during 
realization of the task in question. The proposed solutions 
are mainly addressed to enterprises ofthe above mentioned 
kind. In such enterprises permanent dimensional control is 
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essential, that makes role of measurement processes more 
meaningful and in consequence the discussed problems 
more important. 

e The solutions proposed in this paper are tested by 
GEOMETR company which carries out measurement 
operations in the enterprises producing large-size steel 
structures. 
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